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ABS1RACT 

This study was undertaken so as to develop the observation equations 
necessary to utilize lunar laser ranging and Very Long Baseline Interferometry 
(VLBI) measurements for the establishment of a primary control network on 
the moon. The network consists of coordinates of moon points in the seleno- 
detic Cartesian coordinate system, which is fixed to the lunar body, oriented 
along the three principal axes of inertia of the moon, and centered at the 
lunar center of mass. 

The determination of coordinates of points on the moon using earth-based 
observations requires the knowledge of the following dynamic behavior of the 
earth and the moon: the orbital motion of the moon about the barycenter, the 
rotation of the moon on its axis and the motion of the earth about its center 
of mass. In addition, the knowledge of the geocentric positions of the terres- 
trial stations is essentia:. 

Since our knowledge of the parameters related to the above phenomena can 
be improved simultaneously with the determination of coordinates of lunar 
points, the observation equations derived in this study are based on a general 
model in which the unknown parameters included the following: 

(a) The selenodetic Cartesian coordinates. 

(b) The geocentric coordinates of earth stations. 

(c) Parameters of the orientation of the selenodetic coordinate 
system with respect to a fixed celestial system. 

(d) The parameters of the orientation of the "average" terrestrial 
coordinate system with respect to a fixed celestial coordinate 
system. 

(e) The geocentric coordinates of the center of mass of the moon, 
given by a lunar ephemeris. 
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The orientation parameters of both the earth-fixed and the moon-fixed j 

coordinate system were represented in this study by three Eulerian angles ) 

which, along with their time rates, could be obtained by numerically integrating • 

J 

the differential equations of motion of the respective bodies. This resulted in j 

the reduction of the number of parameters (in the adjustment model), which 

are related to the orientation of the two bodies. The general adjustment model j 

developed for the analyses of laser and VLBI observations is based on the 

theory of adjustment computations with matrix algebra. 

The numerical tests performed in this study with simulated as well as real j 

data demonstrated that the numerical integration of the earth's orientation 
angles yields values which are very close to the classically computed angles, j 

and that the initial conditions for the integration are capable of being solved 
in an adjustment process. Also, the results of the numerical experiments 
performed with simulated laser data confirmed the feasibility of the method 
developed in this study. 
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1 INTRODUCTION 
1. 1 General Discussion. 

The science of lunar mapping is an old one whose origin can be 
traced to the beginning of telescopic astronomy, started by Galileo Galilei 
in the seventeenth century. Since the first lunar map was published by 
Galilei around 1610, numerous authors have been engaged in mapping the 
earth's closest celestial object. Traditionally, this field has been part of 
the domain of astronomers. However, things changed during the advent of 
the space age, and scientists of differing disciplines became interested in 
lunar mapping and related subjects. The quality of published maps of the 
moon improved with time, due to the advancement in both theory and 
instrumentation. A detailed history of lunar mapping has been published 
by Kopal in Chapter 15 of his book on the moon ^51]. 

The mapping of an area requires as a first step, the establishment of a 
consistent network of control points, whose coordinates have been obtained 
in a unified and well-defined coordinate system. As opposed to the method 
of establishing geodetic control for earth mapping, past selenodetic con- 
trols have been obtained through the use of earth-based observations. Thus, 
the determination of coordinates of points on the moon through classical 
methods is rigorously tied to the following extraneous parameters: 

(i) the coordinates of the moon's center of mass in a 
geocentric inertial coordinate system 

(ii) the orientation of an earth-fixed coordinate system 
with respect to an inertial coordinate system 
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(iii) the orientation of a defined set of moon-fixed 
coordinate axes with respect to an earth-fixed 
coordinate system 

(iv) the atmospheric refraction model used for reducing 
the earth-based observations. 

Our knowledge of these parameters is constantly improving due to a 
natural loop of events, ir. which improvement in the various theories 
demands better instrumentation and observation reduction procedures, 
while more accurate observations demand improvements in existing theories. 

The classical types of observations used in lunar position determination 
were the heliometer observations and earth-based lunar photography. In 
many instances, the same observations (with or without additional observa- 
tions) were re-evaluated many times as a result of significant improvements 
in theory. It now appears that further refinements in the reduction methods 
for the old observations can neither significantly improve positional 
accuracy of lunar control points, nor contribute to further improvement 
in theory. Hence, over the past few decades, new forms of observations 
have been used to provide selenodetic control, and the application of new 
instrumentation for selenodetic control has been suggested. The main 
objective of this study is to investigate how two of the new observational 
systems can be applied for establishing a primary network of lunar control. 

In the section that follows, a brief review of past methods of determin- 
ing coordinates of lunar features will be presented. The last section of 
this introductory chapter lists the new observational systems suggested, 
the problem areas to be investigated and the approach used in this study. 


2 






1.2 Review of Past Methods. 

The classical method of determining the moon's rotational constants 
(I - the mean inclination of the moon's equator with respect to the mean 
ecliptic and f - the moon's so-called mechanical flattening), its geometri- 
cal figure and the physical libration uses the reduction of heliometer observa- 
tions. The heliometer is a refracting telescope of small aperture which 
measures angular distances between a "first fundamental" point on the 
moon and other lunar features and the moon's limb. The "first fundamental" 

point is Mosting A, a medium-sized crater situated close to the center of 

) 

the moon's disc. The basic observation technique designed by Bessel in 
1839 has been used with certain modifications to obtain the numerous 
heliometer observations available today (39]. 

Those heliometer observations involving the moon's apparent limb 
and Mbsting A have been used by many famous astronomers to simultaneously 
solve for the lunar physical librations and the coordinates of Mosting A. The 
adjustment procedure involves the fit of a best fitting sphere to the observed 
moon’s limb, and the coordinates of Mbsting A obtained in this manner 
are referenced to the center of figure. 

The determination of coordinates of other fundamental control points 
on the moon utilizes the heliometer observations made between lunar 
features and Mosting A. Using the available physical libration parameters, 
these positions are determined relative to the "first fundamental" point. In 
1899, Franz, a German astronomer determined the coordinates of eight 
fundamental control points on the moon f27l. Later, coordinates of four addi- 





tional craters were obtained by Hayn. Most of the existing coordinates of 
lunar features today are based on the coordinates of these fundamental 
control points. 

The extension of lunar control was achieved through the use of 
earth-based photography. In general there are two types of lunar photography, 
namely phase and full-moon photography, each possessing its own advantages 
and disadvantages. Furthermore, moon photography is either with or with- 
out star background. The star background provides the orientation of the 
camera system as well as the focal length of the telescope. On the other 
hand, photographs without star background rely on the fundamental points 
for scale and orientation. 

With the use of full-moon photographs taken at Lick Observatory at 
different librations, Franz extended the original eight points of Franz into 
a system of 150 points [28l. For the first half of this century this set of 
points was used as the primary control network on the moon. Later, in 
1956-58, an Austrian astronomer Schrutka-Rechtenstamm undertook an 
extensive review of the moon lib ration theory, and then recomputed Franz's 
150 points r 89]. Thus, the Schrutka-Rechtenstamm system was established, 
and has since been used for further control densification. 

In the late 1950's and early 1960's, there was a great need for more 
accurate control extensions on the moon so as to provide a basis for lunar 
cartography and to meet the needs of the manned lunar explorations. Several 
individuals and agencies became involved in the lunar control and lunar 
mapping work, among whom were two U. S. agencies - the Aeronautical 
Chart and Information Center (ACIC) and the Army Map Service (AMS, now 
TOPOCOM). Using different types of moon photography as well as reduction 
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procedures, each agency obtained a set of coordinates for lunar features. 

In spite of large differences between the two systems, a combined 
solution was performed by the Army Map Service in which both the ACIC 
data and the AMS data were used to obtain selenodetic coordinates of 
734 craters. This solution is now referred to as the Department of 
Defense (DOD) Selenodetic Control System 1966 [ 7 ]. 

Another direct result of man's lunar exploration over the past few 
decades is the addition of satellite-borne methods of determining coordinates 
of features on the moon to the classical methods mentioned above. These 
methods include photography of the moon by lunar satellites and direct 
angular observations of features of the moon by astronauts in a manned 
spacecraft. These and other methods are discussed in F 78]. 

The main features of the classical methods of obtaining selenodetic features 
as reviewed above are that they are based on earth-based photographs and 
on a fundamental control network established through heliometer measurements. 
The earth-based photographs suffer from scale and resolution limitations, 
and the heliometer Is a short-focal length telescope with limited resolution. 
Coordinates of points on the moon tied to fundamental control points are 
affected by the center of figure - center of mass bias. Relative accuracies 
of positions of lunar features from earth-based photography are approximately 
1 and 1.5 km near the origin, in the horizontal and vertical components 
respectively, and degenerate rapidly towards the lunar limb f 753. 



Future Observational Systems for Selenodetic Control. 


In the summer of 1965, the National Aeronautics and Space Admini- 
stration (NASA) sponsored a Lunar Exploration and Science Conference in 
Falmouth, Massachusetts [74]. Two years later, a follow-up conference 
was held at Santa Cruz, California [75]. Among the groups participating 
in each of these lunar conferences was a Geodesy and Cartography working 
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group. This group, after clarifying its scientific objectives, went on to 
suggest ways of achieving these goals. The scientific objectives agreed 
upon at the meetings included the following: 

(a) more accurate determination of the position and 
orientation of the moon (and their variations 
with time) with respect to the earth 

(b) determination of the mass and gravitational field 
of the moon 

(c) determination of the physical size and shape of 
the moon, measurements of topographic variations 
and development of a lunar geodetic network over 
the entire lunar surface. 

In order to achieve these goals, the group suggested the use of observa- 
tional systems (some new and others old) such as 

(i) laser ranging (to the moon and lunar satellites) 

(ii) earth-based radio tracking of the moon and lunar satellites (with 
Doppler and ranging transponders on the moon and on the lunar satellites) 

(iii) orbiting camera systems 

(iv) independent-clock radio interferometry (VLBI) 

(v) surface gravimetry, and gradiometer measurements 
taken from lunar satellites. 

To this list was later added satellite altimetry, satellite-to-satellite tracking 
and earth -moon, moon-moon Very Long Baseline Interferometry (VLBI) Tl8]. 

Since these new observational systems were proposed, few papers 
have been presented at different scientific meetings on what some of 
the above instrumentation (notably the laser and VLBI) could achieve in the 
field of selenodesy. However, these papers are usually of a general nature, 
without specific and detailed mathematical development of how these tnstru- 
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ments can be used to improve selenodetic control and related parameters. 

The primary aim of this study is to develop observation equations 
which could be used with new observations to achieve some of the goals 
listed above. In doing this, only the lunar laser ranging and the VLBI 
have been considered. These observational systems offer a 
means by which a new primary lunar control network can be 
established, while at the same time the position and orientation 
parameters of the moon can be Improved. Such a primary control 
network established through the instrumentation considered in this work can 
be extended through other modern methods (see [78]) together with lunar 
surface measurements. The current status of the position and orientation 
parameters of the moon will be discussed briefly below. 

Lunar Position. 

The problem of finding the position of the center of mass of the moon in 
a g.-.ocentric inertial coordinate system has attracted the attention of numer- 
ous scientists since the time of Newton. This has resulted in a few lunar 
theories., each one being an attempt to analytically solve the differential 
equations of motion of a perturbed two-body system. A general account 
of the different methods used in obtaining the moon's motion in space is 
contained in a book by Brown [11 ]. 

The result of any lunar theory is the lunar ephemeri3 which contains 
the positions of the moon in tabulated form. The ephemeris in general use 
today is based on the Hill-Brown theory, with corrections by Eckert C 20]. 
Over the past few years, the Jet Propulsion Laboratory (JPL) has intro- 
duced another method of obtaining lunar position. This involves the numer- 
ical integration of the equations of motion using digital computers. The 



V 

S3S5I (RSSSSSSJCSS? ■ 


7 





i 

i 

) 

I last of the series of numerical lunar ephemeris published by JPL and 

j available on magnetic tape (DE-69), is the LE-16 [ 77] which is believed 

j to be gravitationally consistent and the most accurate ephemeris available 

\ today. Nevertheless the estimated uncertainty in the ephemeris is about 

? 50 meters in range and 100-150 meters in position. The expected pre- 

; cision of laser and VLBI measurements is much higher than that of any 

lunar ephemeris available, and thus the expectation that the new observa- 
i. tions can be used to improve the current status of the limar theory. 

t Moon's Orientation. 

T 

i. Traditionally, the orientation of the moon in space is defined by three 

c Eulerian angles between a moon-fixed set of axes and the mean ecliptic 

£ 

j- coordinate system. The orientation angles are obtained through the solution 

of the moon's rotational equations of motion. The rotation of the moon 
i follows approximately the empirical Cassini laws, and the deviations of the 

\ true motion from that defined by the Cassini laws are the physical librations 

t 

jj of the moon. 

The physical libration of the moon is traditionally solved by applying 
heliometer observations to the linearized forms of the equations of motion. 
Hence the accuracy of these angles depends on the accuracy of the heliometer 
observations as well as the errors introduced in the process of obtaining the 
linearized form of the equations of motion. Currently, lunar orientation 
(i. e. physical libration) is estimated to be accurate to about 20'' of arc 
(about 200-m on lunar surface) [75], 

It has been recently proposed that the lunar orientation angles be obtained 
by direct integration of the original equations of motion T 78], The resulting 
angles and their rates will depend only on the six initial conditions and the 
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moon's physical parameters (for example the principal moments of inertia). 
These parameters can then be corrected in an adjustment using the new 
types of observations. 

Earth's Orientation . 

Since the laser and the earth-moon VLBI are observations involving 
earth stations, the orientation of the earth in space must also be known. 

The earth's orientation in the past has been calculated through precession, 
nutation, and polar motion and its accuracy is considered to be within that 
demanded by current types of astronomical observations. 

The use of modern instrumentation of higher accuracy such as the 
VLBI and laser introduces the possibility of updating our present knowledge 
about precession nutation and polar motion. In order to do this without having 
to adjust all the numerous coefficients of the series expressions for nutation and 
precession, it has been proposed in this work that the orientation of the 
earth with respect to an inertial system be defined by three Eulerian angles. 
These angles and their time rates can be obtained Dy direct numerical 
integration of earth's equations of motion in its original form. The advantage 
of this method is that the earth's orientation parameters are reduced to the 
six initial conditions only. 

In this study, both the earth and the moon have been assumed to be rigid 
bodies, and the equations developed are valid for a rigid earth and moon. The 
effects of non-rigidity of both celestial bodies on the equations derived are two- 
fold. Firstly, there is the direct effect of earth tides and lunar tides on the 
geocentric and selenodetic coordinates of earth and moon stations respectively. 
Some discussion on this effect has been included in Appendix B. The non- rigid- 
ity of the earth and the moon also affects the rotational equations of motion of 
the two celestial bodies in a rather complicated manner. Investigations on 
these two effects arc continuing. 
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The presentation of this work is organized as follows: The next chapter 
presents the equations related to laser measured distances. Partial 
derivatives necessary for the evaluation of observation equations 
are fully derived. Most of the equations in this work are expressed in 
compact forms using either vector or matrix notations. In Chapter 3, equations 
for the VLBI measurements are derived for three cases (earth-moon baseline, 
moon-moon baseline, earth-earth baseline). The numerical integration method 
of calculating the earth's orientation angles is presented in Chapter 4. The 
equations of motion are derived for the case of a rigid earth without any 
other assumptions, and preliminary computational results are included. 

In the fifth chapter, equations for a least-squares adjustment pro- 
cedure for the laser and VLBI observations are given, and preliminary 
computational results using simulated laser observations are presented. 

The last chapter contains the summary of the work done on this study and 
the conclusions reached. 
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2. APPLICATION OF LASER RANGING TO SELENODETIC CONTROL 
2. 1 Historical Background 

Laser, an acronym for light amplification by stimulated emission of 
radiation, is the result of an extension of the principles of microwave ampli- 
fiers (MASERS) to amplification and generation of electromagnetic radiation 
in or near the optical region. With the laser, it became possible to generate 
and control coherent light through the use of electronic transitions in atoms. 

The first successful construction of light amplifier (laser) was made by 
Maiman in 19G0 ^56:. Since then, numerous modifications and improvements 
have been made to the laser, and at present, laser has become one of the 
most important sources or generators of radiation. 

The main properties of the laser are directionality, high intensity, coher- , 
ency.and narrow frequency width (monochromaticity). Radiation generated by 
laser can either be pulsed or continuous. Because of its properties, the laser 
has many applications in numerous scientific fields such as in distance measure- 
ments, communication, medical and biological research, and in the general 
sciences l60]. 

In geodetic science, the laser has been used for the construction of precise 
electromagnetic distance measuring equipment. Also, laser ranging systems 
have been built, for the purpose of satellite tracking and obtaining accurate 
ranges to the moon. Today, a few satellite tracking stations have been 
equipped with lasers, and there is a lunar laser ranging station in the United 
States, actively engaged in ranging to the corner reflectors on the moon. 

In 1965, Alley at al. [ 1 ] proposed the use of a Q-switched ruby laser in 
conjunction with a corner reflector located on the moon for making accurate range 
measurements between the eart'i and the moon. This proposal was independently 
advocated by a group of Russian scientists in the same year [48], A group of 




scientists and engineers was subsequently formed with the purpose of investi- 
gating and planning the lunar laser ranging experiment (LURE). In 1967, the 
National Aeronautics and Space Administration officially endorsed the emplace- 
ment of a retroreflector on the moon as part of the surface experiments to be 
performed on some of the Apollo lunar landing missions. The. first retro- 
reflector was placed on the moon two years later by the crew of the first land- 
ing mission, and the second retroreflector was deposited on the moon by the 
crew of the Apollo 14 mission. A third retroreflector is expected to be placed 
on the moon by the Apollo 15 astronauts. 

Shortly after the first laser ranging retroreflector (LRRR) was placed on the 
moon, reflected laser pulses from the LRRR were acquired with the 120-inch 
telescope at Lick Observatory of the University of California, located at Mt. 
Hamilton, California [2G ]. Subsequently, successful acquisition of reflected 
laser signals were also reported by two other observatories. These are the 
Lunar Laser Observatory of the Air Force Cambridge Research Laboratories 
(AFCRL), located in the Catalina Mountains near Tucson, Arizona, and the 
McDonald Observatory of the University of Texas at Fort Davis, Texas. After 
the initial signal acquisition phase, only the 107-inch telescope of the McDonald 
Observatory is now being used for continuing measurements to both reirore- 
flectors on the moon. So far, a range precision of ±0. 3m has been achieved 
C 4 ], and improvement to ±0. 15m or better is expected shortly. 

The life expectancy of the retroreflector is in excess of ten years T26l. 
Long-term precise range measurements to retroreflectors on the moon could 
not only be used to establish a primary control network on the moon, but could also 
be used for improving our present knowledge of astronomical, geophysical and 
lunar parameters. In fall, Kaula summarized the various geophysical problems 
to whose solution lunar laser ranging might contribute. These problems include 
short- and long-term wobble (of the earth's pole) and variations in the earth's 
rotation, earth tides, tidal dissipation and plate tectonics. However, for the 
successful application of laser ranging to these and other problems, stations 
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on the moon need to be established, while more rctroreflcctors should be 
placed on the moon. 

2 . 2 Basic Operating Principles of 1 aser Hanging Systems. 

The general theory of the generation and amplification of light by stim- 
ulated emission is a subject that is well treated in books on lasers and in 
modern physics textbooks. It is not intended here to give a detailed descrip- 
tion of the laser theory , r.or a complete technical description of laser rang- 
ing systems. Rather, in this section, the basic characteristics of any 
laser ranging system will be treated briefly for completeness, and in order 
to make the equations developed in the rest of this chapter clearer. 

The main components of a laser ranging system consists of the laser, 
transmitting and receiving telescopes and a correlator. 

The laser component serves as a light source, and generates a laser beam 
which may be either continuous or pulsed. The pulsed beam is used in the 
laser ranging systems. Also, all presentlunar ranging stations are equipped 
with Q-switched ruby lasers consisting of an oscillator and three amplifiers. 

The laser pulses are transmitted through a telescope which tracks the 
corner reflector on the moon's surface. A small part of the reflected light 
is collected by the same telescope (or another telescope serving as a receiv- 
ing telescope). The same telescope can be used as a transmitter and a 
receiver partly because the light travel time of about 2.5 seconds between 
transmitting and receiving allows the mechanical insertion of a mirror 
whichdirects the returning photons collected by the telescope into a photo- 
multiplier detector. 

The correlator of the ranging system measures the light travel time 
to the reflector and back. As the laser pulse is fired, a start signal is 
generated to start the travel time counter, and detection of the return 
beam generates a stop signal to stop the travel time counter. With current 
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laser and timing techniques, an uncertainty of ±1 nsec in the absolute 
measurements of the round trip travel time is possible. This is 

equivalent to ±15om in the one way distance. 

Other minor features in some of the laser ranging systems include an 
optical offset guider. This feature is used for pointing to the retroreflector 
when the retroreflector is on the dark portion of the moon, by optically off- 
setting the main telescope from a given identifiable lunar feature. When the 
corner reflector is on the illuminated portion of the moon, the signal-to- noise 
ratio of the return laser pulse is considerably lower for smaller telescopes due 
to the background light gathered from the surface of the moon. 

The laser ranging retroreflectors were designed to serve as reference 
points on the lunar surface to be used for point-to-point distance measurements. 
The two retroreflectors already deposited on the moon by (J.S. astronauts are 
each 46cm-square arrays, consisting of 100 fused silica corner reflectors 
(3.3cm in diameter), each recessed by 1,9cm in an aluminium panel. The 
corner reflectors are uncoated and use total internal reflection. Each corner 
reflector has the unique property that light shining into the corner reflector 
will be sequentially reflected from its three faces and will come out in a path 
parallel to the incident light. The third retroreflector array which will be 
placed on the moon by the crew of the Apollo 15 mission is similar to the pre- 
vious two arrays, but consists of three times as many corner reflectors. 

The design of the LRRR is such that the reflectors will perform under 
essentially isothermal conditions throughout lunar nights and most of lunar 
days, thereby ensuring adequate return of signals at all times. The array 
also allows for the velocity- aberration displacement of up to 1. 6km without 
any significant loss of signal. The most important function of the array, how- 
ever, is that it eliminates the stretching in time of the return signal which 
would otherwise be produced by the curvature and irregularity of the moon's 
surface. 




2.3 Earth-Moon Distance Equations. 


The pri.naiy purpose of this chapter is to investigate how lunar laser 
ranges can be used to establish a selenodetie control network on tne moon. 

The retroreflectors would serve as the "benchmarks", whose coordinates 
could be determined using the lunar ranges from laser stations on t K e 
earth. However, the relative position of a lunar point such as the retrore- 
flector with respect to an earth station depends on a number of physical and 
geometric parameters of the earth, moon and the earth-moon dynamic 
system. It is therefore logical to expect a refinement of these parameters 
in any adjustment of lunar laser distances for the purpose of establishing a 
selenodetie control system. 

Before one can use lunar laser distances for the purpose mentioned 
above, the mathematical relationship between the measured distances and 
the following list of parameters have to be established: 

(a) The selenodetie coordinates of the lunar retroreflector (M) 

(b) The geocentric coordinates of the selenocenter (C) 

(c) The geodetic coordinates of the laser station (P) 

(d) Parameters of the orientation of a moon-fixed coordinate 
system with respect to an inertial coordinate system. These 
are given by the physical libration angles < 7 . p. r and the mean 
orientation angles as defined by the Cassini laws 

(e) Parameters defining the orientation of the "average" terrestrial 
coordinate system with respect to an inertial system. In 
Classical Astronomy, these parameters are usually given by 
the processional elements (£o. 9 if ?.), the nutational parameters 
(Ae, £$>), time (represented by Greenwich Apparent Sidereal 
Time (GAST)), and the position of the instantaneous true pole 
with respect to the Conventional Internation Origin (CIO pole), 
given by x p and y f fG9]. 



From the review of existing literature on this subject, one could find 
that a few mathematical formulas have been derived, which express the 
earth-moon distance as a function of some of the parameters listed above. 

For example. Alley and Bender treated the problem of using laser distances 
to improve our knowledge of geocentric longitude in T 2 ], while Alley et al. 
offered a simplified model for estimating the expected accuracies with which 
some geophysical and lunar parameters can be determined in [3 ]. A 
relatively more detailed mathematical model was developed by Kokurin et al. 
for improving certain astronomical parameters in [48] and [49]. Their 
approach was based on geometric consideration of the earth's and lunar 
positions at the epoch of observation. In this section, a different approach 
will be used, based on the consideration of the various coordinate systems 
involved, which can be transformed to a uniform coordinate system through 
the use of rotation matrices (composed of functions of the orientation parameters) 
and translations. The main goal will be to express the coordinates of the lunar 
point at the observation epoch, in a coordinate system topocentered at the 
observing station. The earth-moon distance can then be deduced as a function 
of these topocentric coordinates. 

Before going into the development of the laser equations, the various 
coordinate systems involved and the parameters connected with them will be 
given brief treatments. 

2.31 The "Average'' Terrestrial Coordinate System . 

The average terrestrial coordinate system is defined by the average ter- 
restrial pole of 1900-05 (designated as the Conventional International Origin 
(CIO)), and an average terrestrial equator r 69]. It is the system to which 
absolute geodetic and reduced astronomic coordinates of any physical point on 
the earth are referred. 
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The average terrestrial Cartesian coordinate system (U,V,W) has its 
origin at the center of mass of the earth, and its axes are fixed with respect 
to the earth's body. The positive W axis is directed toward the CIO pole, 
the positive U axis is directed towards the Greenwich Mean Astronomic 
meridian as defined by the Bureau International de l'Heure (BIH), and the 
V axis completes a right-handed rectangular system. 

In computing geodetic coordinates of points on the earth's surface, a 
reference figure (the reference ellipsoid) is usually introduced. An absolute 
geodetic coordinate system can be defined, which is geocentered and has its 
axes coincident with the U, V,VV axes defined above. If the absolute geodetic 
coordinates of an earth station is given as (tp, A, h), then the Cartesian 
coordinates can be computed from: 


‘u~ 


(N + h) eos<p cosX 

V 

= 

(N+h)cos<p sin A 

w 


(N(l-e^ + h) sincp 


where 

N = 3 . 2 -: i , the transverse radius of curvature, 

(I - e sin <p)f 

e a = 21-f 3 , the first eccentricity, 

and 

a, f are parameters of the earth-ellipsoid, denoting the semi- 
major axis and flattening of the ellipsoid respectively. 

2. 32 The Selenodetic Coordinate System . 

On the moon, we use spherical coordinates - longitude -t, latitude b and 
radius r - to define the location of a point on the lunar surface. Similar to 
the geodetic coordinate system, the selenodetic coordinate system is fixed 
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with respect to the lunar body. 

The selenodetic latitude is measured from the lunar equator, defined as 
the plane passing through the center of mass, perpendicular to the axis of 
rotation. Latitude is measured in degrees positive toward the north lunar 
pole and negative toward the south pole. 

The prime meridian, -t =0 is a plane which contains the axis of rotation 
and the Earth-Moon line when the node and the perigee coincide (i. e. at 
zero geometric libration). The longitude is measured in degrees, positive 
in the direction of rotation (eastwards) from 0° to 180° and westwards from 
0° to -180°. 


z North 



Figure 2. 1 The Selenodetic Coordinate System. 
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The Cartesian selenodelic coordinates of any lunar point (denoted *-v 
x,y ./.) is related to the spherical coordinates of that point by the liquation: 


| N '! 

r 

cos b 

cos f. 

>■ 1 - 

r 

cos b 

sin L 


r 

sin b 


The xy 7 . system is centered at the moon’s center of mass. The /. axis is 
coincident with the moon's axis of rotation, the x axis coincides with the 
mean Karth-Moon line (sometimes referred to as "first radius") and ’he 
y axis is perpendicular to both x and z axes and completes a right-handed 
rectangular system. 

The sclenodetic coordinate system is related to the ecliptic system by 
the three tvulerian angles, which will be treated in more detail in section 2. "4, 

2. 33 The Geocentric Coordinates of the Selcnoccntcr. 

In developing the equations for measured lunar laser distances, the 
coordinates of the lunar center of mass with respect to the gcocentcr arc 
needed. These coo xiin.ites are obtainable from any lunar ephemeris which 
gives the coordinates of the moon, relative to the earth, as a function of time. 
All ephemerides in present use are the result of a particular lunar theory which 
attempts to apply the law of gravity to the motion of the moon. 

The lunar theory is a subject that has been discussed and treated by var- 
ious scientists from the time of Newton to the present day. Although different 
authors employ different methods in arriving at their lunar theory', they 
all start from the differential equations of a perturbed two-bodv system. 

The main problem is to integrate the equations of motion of a perturbed 
two-bodv system given by 

f v G (Alp *■ M„)-p = 
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r is the earth-moon vector 
G is the gravitational constant 

M®, Mj are the earth's and moon's mass respectively 

and 

R is the "disturbing function". If there were no 
perturbing bodies, R would he zero and the result- 
ing motion would be elliptic. If the sun acts as 
the only perturbing body, then R is given by 



where 

pis the moon-sun distance 
r' is the earth-sun distance 
Mq is the mass of the sun. 

The disturbing function generally used in any of the lunar theories is modified 
to account for the presence of other perturbing bodies, i.e., the planets. 

The analytical integration of the equations of motion for thn i or more 
bodies in a closed form is not generally obtainable. Consequently, there are 
two basic ways in which the lunar ephemeris is obtained: 

In the first method, the elements of the equations of motion are expanded 
as Fourier series, yielding a sequence of terms that are analytically integrated, 
and resulting in lengthy analytical expressions for the lunar ephemeris. In 
the alternate method the original unexpanded equations of nr ition are numeri- 
cally integrated. The analytical method has been used by classical lunar 
theorists and its main disadvantage is that the expression for the lunar ephemeris 
is truncated. The number of terms in the final expression depends on the number 
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of bodies and types of perturbations treated as well as on the accuracy of the 
high-order terms required. The present most widely used lunar ephemeris - 
which is the result of the Hill-Brown theory [ 12] with modifications by 
Eckert T20] - uses the Fourier series development method. 

The numerically integrated lunar ephemeris, which was initiated by the 
Jet Propulsion Laboratory f.TPL), is made possible by the advent of digital 
computers and is noted for ita gravitational consistency. All the planetary 
effects are well taken care of and any differences between the ephemeris and 
the true coordinates of the moon are a direct result of inadequacies in the 
theory itself. The numerical ephemeris is also clearly superior to the 
analytically obtained ephemeris in terms of accuracy as demonstrated by 
the JPL [ 72]. Consequently, the position of the center of the moon, to 
which corrections are to be computed and which should be used in the pre- 
diction equations for laser distances, should be that from the latest of the 
numerical lunar ephemeris - LE16 which is contained in the JPL Develop- 
ment Ephemeris Tape 69 (DE-69) [ 77], 

2. 34 Parameters of the Orientation of the Moon. 

The laser distances, like any other earth-moon observation, is a function 
of the parameters that define the orientation of the moon. The orientation 
of the moon is related to the moon's rotation on its axis in the same way as 
the lunar ephemeris is related to the orbital motion. 

The moon's rotation follows, to a high degree of approximation, the 
empirical Cassini larvs. The physical deviation of the moon's motion 
from the described by the Cassini laws is known as the physical libration 
of the moon. 

The moon' 8 orientation in space is defined by the three Eulerian angles, 
which define the orientation of a moon-fixed coordinate system with respect 
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to an inertial system. Traditionally, the inertial system chosen is the 
mean ecliptic coordinate system of date which is not truly inertial, but 
has secular changes due to planetary precession. 
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Figure 2.2. The Moon's Eulerian Angles. 


Figure 2. 2 shows the xyz coordinate system which is fixed to the body of 
the moon, and the XYZ system representing the mean ecliptic system of 
date. The Eulerian angles are denoted by 8, $ and 4. 

The Eulerian angles are composed of two parts, namely the mean angles 
(which are the mathematical expressions of the Cassini laws) and the pertur- 
bations (given by the physical llbration angles). The angles are given by 
the following expressions: 

$> = L+ir-0 + r - a 
t!) = O + a 

8 = 1 + o 
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and 


L is the mean longitude of the moon in its orbit 
Q is the longitude of the mean ascending node of the moon’s orbit 
I is the mean inclination of the lunar equator to the ecliptic 
T is the physical tibration <a longitude 
a is the physical libration in node 

p is the physical libration in inclination. 


The mean angles L and Q are obtainable as explicit functions of time from 
the following expressions taken from [65]. 

L = 270?434 1639 + 481,267? 88.3 1 117 T 

-0? 113 3333 « 10 a T a + 0?188 8889 * Iff 6 T 3 

0 = 259? 183 2750 - 1,934? 142 0083T 

+ 0? 207 7778 * 10*1* + 0? 222 2222* 10 s T 3 


where T is the time in Julian centuries of 36525. 0 ephemerls days from Jan- 
uary 0.5 1900 {J. D. = 2415020. 0). The mean inclination I is given by; 


I = 5521 ?5. 

The physical libration angles are expressed as harmonic series having 
fixed coefficients and arguments which are expressed as functions of time: 

r = C A, cos (a, t + b.) 


with similar expressions for a and p. 

The Eulerian angles which define the orientation of a set of moon-fixed 
axes with respect to an inertial system are obtained through the solution of 
the equations of motion for the rotation of the moon. Classically, the dif- 
ferential equations (represented by the Euler's dynamical equations) are 
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analytically solved by series expansion, and actual observations (usually 
heliometric observations) are used to solve for the coefficients of the harmonic 
expressions. Recently, the numerical integration of the differential equa- 
tions of motion in their original form has been proposed, [78], Byusingthis 
method, the Eulerian angles would be obtained directly by integration without 
splitting the angles into two parts - the mean angles and the physical libration 
angles. Inaccuracies in the orientation angles obtained by numerical integration 
can only result from inadequacies in the lunar rotational theory, and not in 
the type of linearization and approximations that are essential for solving 
analytically the differential equations of motion. 

2. 35 The Parameters of Earth's Orientation . 

In observations made from an earth point to a moon point, the relative 
orientation of an earth-fixed coordinate system with respect to a moon- 
fixed coordinate system is needed. The coordinates of the earth station 
are expressed in i n earth-fixed system while the moon station coordinates 
are in a moon-fixed system. In practice, the problem of determining the 
relative orientation of the moon with respect to the earth at any epoch is 
solved by finding the orientations of the earth and the moon with respect to 
an inertial system. 

The orientation of the earth in an inertial coordinate system is also needed 
in geodetic astronomy, where observations are made to celestial objects, whose 
positions are given in a celestial coordinate system of a standard epoch. The 

precessional elements (To 9i and z) define the orientation of the mean celestial 
equator and equinox of date with respect to the mean equator and equinox of an ini- 
tial epoch. The orientation of the true equatorial coordinate system to the mean 
equatorial system is also defined by the nutation in longitude and obliquity (Ai/i, Lt). 
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Lastly, the orientation of the average terrestrial coordinate system with respect 
to the true celestial equatorial system is given by the two coordinates of the 
true pole from the average (CIO) pole, together with the Greenwich hour angle 
of the true vernal equinox (Greenwich Apparent Sidereal Time). 

The expressions for precession and nutation elements are obtained from 
the analytical solution of the differential equations of motion representing the 
rotation of the earth on its axis. The nutation series that are currently in 
use in astronomical work are those obtained by Woolard 194.1. The accuracy 
of the precessional and nutationai elements obtained through the use of 
current expressions are consistent with the accuracies of present day observa- 
tional systems. However, in using observations made by future systems 
which are expected to be more accurate, there is a possibility of improving 
the accuracies of the elements of precession and nutation. The constants 
for the expressions that give Co . 9i and Z; as veil as those for the nutation 
series can be regarded as pai ameters to be adjusted in addition to other 
parameters of interest. 

Another way in which the ortentation of an earth-fixed coordinate system 
can be obtained from the earth's rotational theory is proposed in this 
work. Essentially, this method considers the orientation of the terrestrial 
coordinate system with respect to a fixed mean ecliptic system as made up of three 
Eulerian angles. The three angles and their time derivatives can be obtained 
from the three second order differential equations of motion for the rotation 
of the earth in its complete form by numerical integration process. This 
topic is being given a complete treatment in Chapter 4. 

Figure 2. 3 shows the relationship between the average terrestrial 
coordinate system (UVW) and a mean ecliptic coordir.ate system (XYZ) of a 
standard epoch. The earth's Eulerian angles are denoted by 8 f , >j) t , and 4> t , 
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Figure 2. 3 The Earth's Orientation Angles. 


The advantage of the numerical integration method of obtaining the 
orientation of the earth over the classical analytic method is the fact that 
the angles obtained depend only on the initial conditions represented by six 
constants. Consequently, only six parameters need be added to the param- 
eters list if the values of precession and nutation elements obtained from 
existing expressions are to be taken as variable quantities. 

2. 36 Computation of Distances between Points on the Earth and on the Moon. 

In Figure 2.4, P is a station on the earth miking distance measurements 
to a moon point M. The geodetic coordinates of P are given by <?, X and h. 
The Cartesian coordinates of P in the UVW system are given by 
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Similarly, we let l, b, r be the selenodetic coordinates of point M. Then 
the Cartesian coordinates of M in the xyz system are given by 

I x J] fr cos b cos -O 


y„ = r cos b sin l . ( 2 . 2 ) 

z„ r sin b 

The geocentric Cartesian coordinates of the lunar center Xc , Y c , and Zc at 
the epoch of observation is taken from any lunar ephemeris. In the deriva- 
tions to be made here, it will be assumed that the Xcq, Ye,. Zc q coordinates 
of the selenocenter is obtained from the latest numerically integrated lunar 
ephemeris - the LE16. In this ephemeris, the coordinates are expressed in the 
1950. 0 mean equatorial system. 

The coordinates of the laser station P can be expressed in the 1950. 0 
mean ecliptic coordinate system by applying a series of orthogonal trans- 
formations [ 69]: 


”xTj 

Y 'l 

_ z Jrcu 


= RjUoJp'nV 


where P, N, S are 3*3 orthogonal transformation matrices as follows 
P = Ra(-z0 Raft) RaKo) 

N = Ri (-€ - Ac) R 3 (-A$) Rj (e) (2.4) 

S = R s (-Xp) Ri(-y P ) R 3 (GAST). 

The conventional orthogonal rotation matrices R,(or), R 2 (cr) and Ro(a) are 
often used (as in this case) to rotate a general system by angle a about the 
axes X,Y,Z respectively. These matrices, consistent with a right handed 
coordinate system are given by: 
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10 0 
0 cos a sin a 

0 -sin a cos rv 


R 2 (*) 


cos a 0 -sin a 
0 10 
sin a 0 cos a 


R 3 (») 


cos fy sin a 0 
-sin a cos a 0 
0 0 1 


The parameters z t , 8j and £o are the precession elements for the period between 
1950.0 and the epoch of observation; Ac and A/) are nutation in obliquity and 
longitude respectively at the epoch of observation; e is the mean obliquity of the 
equator of date; CAST represents the Greenwich Apparent Sidereal Time and 
x p . y p are the polar motion parameters. c 0 Is the mean obliquity of the equator 
at the epoch 1950. 0. Equation 2. 3 could be equivalently written as 


'Xp' 


u’ 

Yp 

= Ri(c 0 )Ra(Co)R2ea)R3(z)Rie^R3(Ai»Ri(€ + Ac) Rs(-G AST) R x (y*) Rs(xp) 

V 

Zp 

(Cl. I960. O 

w_ 


(2.5) 


In a similar manner, the coordinates of lunar retroreflector (M) can 
be expressed in the 1950. 0 mean ecliptic system as follows: 
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In equation (2. 6), $, 9 and <J> are the Eulerian angles which define the 
orientation of the selenodetic (x,y,z) coordinate system with respect to the 
mean ecliptic system of date. The three angles have been treated in Section 
2. 34. 

If 

Ti = Ra(4) Ri(-9) Ra(ii>) 


then equation (2. 6) becomes 


XT 


Xm 

Y* 

= R 1 (e e )P T Rx(-c)Tl 

y» 

z„ 

CC1. 1950.0 

_Z* 


(2.7) 


Since the geocentric mean equatorial (of 1950. 0) coordinates of the lunar 
center as obtained from the LE-16 or any other lunar ephemeris are given 

by 

~*r 

Y;q , 

_ Zc Q_ 

the geocentric coordinates of the selenocenter, expressed in the mean 
ecliptic 1950. 0 system are given by 


1— "”1 


XCQ 

Y e 

= Ri(fo) 

Ye, 

_Zc. 

Eel. 1950.0 

_Zcq_ 


where c 0 as defined earlier, is the mean obliquity of the equator at 1950. 0. 

Thus, the topocentric coordinates of the lunar reflector (M), with the 
laser station (P) at center are obtained, by considering Figure 2.5, as 
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topoe. Eel. 1990.0 - J |el. 1960.0 


Eel. 1960.1. 


The topocentrie coordinates of the moon point obtained by equation (2. 9) are 
expressed in the mean ecliptic coordinate system of 1950. 0. It should be 
noted, however, that since distances between two points are independent of 
the coordinate system, the topocentrie coordinates could have been expressed 
in any coordinate system. The 1950. 0 mean ecliptic coordinate system has 
been chosen because it is a near-inertial coordinate system which will also be used 
in the case of VLBI observations involving directions. 



Figure 2. 5. Topocentrie Mean Ecliptic (1950. 0) Coordinates 
of Retro reflector (M). 
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Into equation (2.9), we substitute expressions for 
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which is fully expressed as 
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where 

X - is position vector, the subscripts indicating 

T - moon point in a topocentric coordinate system 
CQ - geocentric equatorial coordinates of moon's center 
of mass 

M - moon point in selenodetic coordinate system 
P - earth observing station in average terrestrial coordinate system 
M lt M a , M 3 are orthogonal transformation matrices for expressing 
each position vector in a coordinate system parallel to that of the 
1950.0 mean ecliptic system. 

The rotation matrices in equation (2.5) were used in order to transform 
from the average terrestrial system to the mean ecliptic system. It has been 
explained in Section 2. 35 that the orientation of the average terrestrial system 
to the 1950. 0 mean ecliptic system could be defined by a set of three Eulerian 
angles - 8 £ , ipt> and . If the three orientation angles are available, then 
equation (2. 5) can be written as 


( 
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Cel. 1950.0 


= T 3 t V = RaH>£) Rife) Rat- *£> V 


T a = i^(4t)Ri(-e t )R3 (<M ( 2 * 13 > 

In a similar manner, the numerically obtained Eulerian angles for the 
moon's orientation (see Section 2. 34) could define the orientation of the 
selenodetic coordinate system to an inertial system such as the mean 
ecliptic system of 1950. 0. Therefore equation (2. 7) could also be written 
as 


y« = RafiW Rife) R3 {“^r) y» 


ECU 1003.0 

where 0 O , 0„, are the numerically integrated moon's Eulerian angles 
defining the oriencation of the selenodetic system with respect to an inertial 
system (in this case, the 1950. 0 mean ecliptic system). 

As an alternate expression to equations (2. 10), equation (2. 9) could be 
expressed as follows: 


~xj 

x " > i 


Y« * = Ri(€o) 

Yt, i 

+ i(3(-</)«) Ri(9r) Rs("4ii) y« 

z„J 

tppoc. 


i 

L Z| L 


- ffcH>t)Ri(8E)R3<-&) v 
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Equations (2. 10) express the topocentric coordinates of a moon point 
at any epoch in terms of astronomical parameters in current usage In most 
astronomical work. On the other hand, equation (2.15) gives the same 
topocentric coordinates of the moon, but expressed in terms of Eulerian 
angles defining the orientation of the earth-fixed and the moon-fixed coor- 
dinate systems with respect to the 1950.0 mean ecliptic coordinate system. 

The distance from the earth point to the moon point is a function of 
the topocentric coordinates of the moon point as follows: 

D = <X£ + Y* + Z%) & (2.16) 

where D is the distance and the X^,, Y m . are topocentric coordinates of 
the lunar point at the epoch of observation given by any of equations (2. 10) 
or equivalently, equation (2. 15). 

2.4 Formulation of Observation Equations for Earth-Moon Plntancpa 

In Section (2. 36) a set of equations have been derived, by which the dis- 
tance between any point on the earth and another point on the moon can be 
predicted for any epoch of observation. In general, the predicted and the 
observed values of the distance will not be identical due to many factors. 
These include the fact that some or all of the parameters used in predicting 
the distance may be inaccurately known, and the fact that the observed 
distances themselves may suffer from systematic errors, in addition to the 
common accidental observational errors. 

In this section, the effect of small variations in the parameters on the 
calculated distance will be derived, which together with the differences in 
observed and calculated values, could be U3ed to correct initial values of 
the parameters. It is assumed here that the measured distances have been 
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corrected for all instrumental and other systematic errors (such as 
atmospheric refraction). Consequently, no systematic error models are 
introduced into the observation equations. Rather, the equations can serve 
as a general model which can be used to adjust any type of distance observa- 
tions from an earth point to a moon point if such observed ranges are freed 
from systematic errors. 

From equation (2. 16), the distance equation can be written as 

D = f 1 ([X t l) (2.17) 

where D is distance and [X, ] is column vector of topocentric coordinates 
of the lunar point. 

Also, from equation (2. 10), if GAST = © 


[X T ] = f a (Co.€.Co.0i.zi. &!>,&(, 6, Xp,y P ,ijb,9,$, U,V,W, 
y*« Z*, Xcq« ^CQ* ^Cq) 


(2. 18) 


and if the parameters listed in equation (2. 18) are represented by a column 
vector [ xl, then 

[X T ] = feta. (2.19) 


The mathematical structure represented by equation (2. 17) with equation 
(2.19) is of the type 

I* = F(X») (2.20) 


where L, is a vector of the adjusted values of observed quantities, and X. 
is the vector of adjusted values of parameters. The Taylor series lineari- 
zation of (2. 20) leads to the observation equations [94] 


where 


V = AX + L 



( 2 . 21 ) 









X = X, - Xo 
L = Lq - Lfc 
Lo= F(Xo) 


( 2 . 22 ) 
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and 

Xo is the vector of initial values of parameters 
Lj, is the vector of observed values 
V is the vector of residuals. 

In this case, the design matrix [A] is given by: 

. _ 3D ap AX, 
ax ax t dx 


(2.23) 


where P. x, X x are all column vectors. For each measured distance d„ 
the ith row of matrix A is given by 


ad, 

OX 


adj aXj 

ax T ax 


l-'rom equations (2. 16) and (2. 17), 


MX T ) = (X^ + + Zg T )^= P. 


Therefore 


and 


^4 = i- [xJ : Y ‘ 

a3? T d, l ^r . x "r 


ap . 
ax T " 


i r < 


i r , 








where n is the number of observations. 
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(2. 24) 


(2. 25) 


(2. 26) 



In order to find — T , it can be seen that the function f a in equation (2. 18) 
includes rotation matrices, which also have to be differentiated. Differ- 
entials of rotation matrices are easily obtained through the use of special 
skew matrices (hereafter referred to as Lucas matrices), whose properties 
are such that when the skew matrix is either pre-multiplied or post-multi- 
plied by a rotation matrix, the resulting matrix is the partial derivative of 
the original rotation matrix [58]. 

Define: 



Then as an example, if R„(a) represents a 3 « 3 rotation matrix, (n = 1,2,3), 


dR,(a) 

-dr * I* - R»<«) = R 0 (O) • Lb 


3Rn (-QQ 


da 


-Lb - R(a) 


(2.27) 
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Also, if 




M = R D (ot) R a (S) Rj(y), 


= L B R B (a)R,f(9)R J (y) = L b M 


^ M 

fj = R,(a)L.R,(fl)R,(y) 




= R B (a)R B tf)R j fv)L J = M Lj . 


(2.28) 


The partials - — t will now be derived from equation (2.10), for each 
o w 

element of the x vector individually: 

(1) For the mean obliquity of the equator at the epoch 1950. 0, (t 0 ): 



/k.i 


*. = LiRxfcc) Ycq + P’Tl y„ - PVS* V 


- Li X T - Lj 


VI 

Y„ 


(2. 29) 


(2) For the mean obliquity of the equator of date (c): 
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(5) For the moon's Eulerian angles, 0, and 


jtf = Ri(fo) p’RieoRne^uRxtej^e^ y * (2.36) 


= -Ki(Cd)P T R } ec) La tI y* 


(2. .17) 


= -RAfo)P T RiH)Tl I* y* . 


(2. 38) 


(6) For the coordinates of the true pole x P , y p and the Greenwich Apparent 
Sidereal Time 0: 


= -RiieoJP'^S 1 !^ V 


(2. 39) 


= -R:(C 0 )P T N T R,(-e)L l R;(y,,)R z (x 1( ) V 


(2.40) 
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(2.41) 


(7) In order to find the partials with respect to the coordinates (of earth 
point, moon point and moon's center), let the column vectors Xeg, X<« and 
Xp be given as (see equation (2. 10)): 


X co = : Y co , X* = y„ . X, = V 


I = identity matrix = 0 1 0 


0 0 1 


1 0 0 


Ri (fo) * l - Ri (<o) • 


(2.42) 


Similarly, 


R 1 (Co)P t Ri('OT 1 t 


(2.43) 


= -R s (fo)P T N T S T . 

oXp 


(2.44) 




I-'-- . 




Each of equations (2.42) to (2.44) is a 3 x 3 matrix, with the three columns 
representing the three coorilinates (of selenocenter, lunar point and earth 
station), while the rows represent the three topocentric coordinates of the 
lunar point. 

In using measured laser distances as observed quantities it will be 
impracticable to use some of the parameters that have so far been used in 
finding the partials represented by equations (2.29) to (2.44) due to one 
main reason. Some parameters, in their present form, change from epoch 
to epoch and the total number of parameters to be solved for will always be 
larger than the number of observations. Such parameters that vary with time 
include the precessional elements, the mean obliquity of the equator of date, 
nutation in obliquity and longitude, the Greenwich Apparent Sidereal Time, 
the Eulerian angles of the moon and the geocentric coordinates of the seleno- 
center. It is therefore necessary to express these variable parameters as 
functions of time or like arguments, and to regard the coefficients of these 
arguments as the parameters. 

Suppose a variable parameter, q can be expressed as a function of time: 

q = go ♦ &t + g 2 t 2 + • • • + g n t n 

i.e. 

q = f(G,t) 

where 

G is a column vector of coefficients go, g x — 



and t is some time unit, or similar argument. 

Then 

ax t = ax, 

bG b q bC 

Each variable parameter will now be considered. 
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! An expression for the mean obliquity of the equator of date (e) is given 

| by T 65]: 

i 

j 

■ € = 23°. 4522944 - 0?01301250t - 0?1638889 x 10~ s t s + 0?5027778 x 10*1? 

( 

| ° r 

j € = 3o + Sit + Hat 8 + % t a 

I , 

where t is in tropical centuries from 1900. 0 

ax T = ax* _ a ( 

oc dla^ • • • ag] 1 

’ and 




a< 

' • * *<3 lT 


[1 


t 2 : t 3 ]. 


If the variable parameter eis then replaced by the coefficients Sq, a r . a^ as. 
then equation (2. 30) is replaced by •. v 


»G^r It ’-' 1 * •* ei 


where — - 1 is a column vector given by equation (2. 30). The right hand side of 
o € 

equation (2.45) is a 3 x 3 matrix which will have to be pre-multiplied by equa- 
tion (2. 25) so as to obtain a row vector that represents the partials of the 
lunar distance with respect to the coefficient Sq, ••• 83 , i.e. , 

ad, = ad. _ dXj _ ac 

oTag, •••a 3 3’ ax| dc a r ao, • • • a ,] 1 * > (2.46) 

, j ) 

The Precessional Elements. 

' VC 

The expressions for the precessional elements are given [ 65]: 








Co = 2304" 951665 1 + 0."302165t s + O^OlSt 3 

6i = 2004?258257 1 - O.^OSSSt 2 - 0."0418t 3 

Z i = 2204 .* '95 1 615 t + l."095195 1® + 0."01832^ 






From equations (2. 47), (2. 48), (2. 49) and (2. 25) the partials of the distance with 
respect to Bi, and Bg can be obtained In a manner similar to equation 
(2.46). ' - ' V. 

The Nutattonal Parameters. 

When considering the nutational parameters Ae and A$, the situation 
is more complicated. The series normally used in calculating Ae and A0 
for each epoch contain some 40 terms for Ac of which 24 terms are of 
short period and 69 terms for A0 including 46 short periodic terms [65], 

To Include all the 109 coefficients as parameters in a general adjustment 
model not confined to solving for only the earth's precessional and nutational 
terms would necessitate a large number of observations as well as increase 
greatly the dimensions of the resulting normal matrix. Moreover, the 
argument in each of the terms of die nutation series is a combination of 
trigonometric functions of the five so-called "mean longitudes" (of the sun's 
perigee, the mean's perigee, sun, moon's node and the moon). The 'mean 
longitudes" are also calculated as functions of time (see [65]). The coefficients 
of the terms in the expressions for the mean angles could therefore be also 
treated as variables. : -X\ V-- ■ 

There are two possible ways in which the problem of the large number 
of parameters can be avoided. In the first alternative, the mean longitudes 
can be regarded as constants, while only a few of the coefficients of the argu- 
ments in the nutation series are included as parameters. The other coeffi- 
cients will have to be held fixed. The choice of the coefficients to be included 
as parameters is more or less arbitrary, but should include coefficients of 
both short and long periodic terms. The second way is to replace the orienta- 
tion of the earth-fixed coordinate system with respect to an inertial system 
by the three Eulerlan angles mentioned earlier, thereby combining the effects 
of 1 uni-solar precession, nutation and planetary precession. In this case, the 
angles can be Integrated numerically and the number of parameters is reduced 




to only the six initial conditions for the three second order equations of motion 
(see Chapter 4). Expressions for both methods will be given here. 

Let the series expressions for Ae and Ati> be given as 

*o 

= S Jj Gt 
1=1 

69 

*ili = £ K, Q. 

i=i 

where 

J,, K, are numerical coefficients 

G t , Q t are arguments (functions of the "mean angles"). 

Also, suppose only n and m coefficients are to be regarded as parameters 
from Ac and A# respectively, such that the reduced expressions are 

A** = £ J, G, 

i=i 

A0* = £ K, Q, . 

1=1 

Then 



= o 7 

iv H 


where 

J, K are column vectors of dimension n and m respectively 
and G, Q are column vectors of dimension n and ir. respectively. 
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Thus 


ax T _ ax T t 
aj aA c 


= Ml 0 t 

aK d&tp ^ 


(2. 50) 


(2.51) 


Equations (2.50) and (2.51) can be used together with equations (2.34), 

(2. 35) and equation (2. 25) to find ‘he partials of the measured distance d t 

4 

with respect to the nutation parameters. The equations for the partials, 
when Eulerian angles are used will be considered later in this chapter. 


The Moon's Eulerian Angles. 

The moon's Eulerian angles 9, ip and 4 are conventionally computed 
through the equations given in Section 2. 34. Here, the mean angles L, O 
and I can be regarded as exact while the libration angles r, a and p are 
variable parameters. Hence, from the expressions for fl, ip and 4 

ax T ax T a4> 

3T a4 ' dr 


ax* _ ax T _ ai|) , ax T a4> 
a a dip ' da d4> 3a 


ax T _ ax., a 9 
a P 30 dp 


(2.52) 


However, since r, a and p are calculated through series expressions, the 
partials of X x with respect to the parameters in the series expressions are 
the required partials. The situation here is similar to that of nutation 
elements. The expressions for T, a and p contain 18, 11 and 11 terms respec- 
tively [65], Some or all of the 40 fixed coefficients can be regarded as param- 
eters if the arguments are regarded as known, or if numerically integrated Eulerian 



EBZaAM 


Hi ’ jwjjarrt.c 




i • .-• --^v- ' 

' -^^ST^ gj^g S ^ g S ^g ESg a aes a EasaE i^BgagBesaa ga^M^gBMMM 


l' 



»u 

P: 


|- 


£ 

e- 

»•' 

I 

i 


v. 

r 


& 

l 

I 

'€_- 

% 

W,’> 

if 

V'V. 
f •/ 1 

c'rv 

at. v- 
Sl'r'J 

&*•' 

fe- 

<«•:,- 

P#‘ 

?#£ 

S3?, 

S*j?7 


angles are used, the six initial values of the integration are the only param- 
eters that need be considered. 

The physical libration in node, inclination and longitude have the respective 
forms: 

1 11 

o = 7 s Sj 

1 1=1 

n 

p = £ F, H, 

t-i 

la 

T = £ G, U. 

i=i 


where E t , F, and G, are constants, S f , H, and U t are sine functions of com- 
binations of the Delauny variables and I is the mean inclination. Thus, if E, 

F, and G represent the column vector of parameters (related to physical libration 
angles ) to be adjusted, then 


d cr _ J_ 
dE ~ I 



11 

dG 


U’ 


( 2 . 53 ) 


where 


S, H, U are the column vector matrices of arguments. 


Furthermore, the mean inclination I can be regarded as a parameter, hence 


!£ 

di 


s 


and from Section 2. 34 



( 2 . 54 ) 


( 2 . 55 ) 
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Thus 






ax, ax, a<t- ar 

dG a$ ar ‘ aG 

ax, . ax T ai/) ag , ax, o^> ag 
aE ai») ' da ' aE a* * a<j * aE 


ax, ax, ae ap 
aF “ afl ‘ ap ‘ aF 

and 

ax, _ aXj ap | ax t aft ag t ax, ag 

ai as ‘ ai a* ’ a * ai a ib ’aa’ ai* 

Noting from the expressions for the moon's Eulerian angles given in Section 
2.34 that 

aft a<& _ ae ae 

ar = da dp ai 

and 




the above equations can be written as: 



ax. 

Mr 

dr 




aG 

aft 

■ d G 


(2.56) 

ax, 

ax T 

da 

ax, _ 

ag 


a e 

dip 

' oE " 

dV 

aE 

(2. 57) 


ax T 

aF 

ax, 

ae 

. & 
aF 


(2.58) 


ax. 

ax. 

ax, ag + ax, _ ag 

(2.59) 

ai 

ae 

" aft ‘ di d</> ' ai 
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ft 

P: 





. v • ■ . ■ * I * 


I 

Equations (2. 56) to (2. 59) can be evaluated because expressions for the j 

partials appearing on the right hand side of the equations have already been j 

given. The size of the matrix obtained from equations (2. 56), (2. 57) and ^ 

(2.58) remain 3 * n, 3 * m and 3*4 where n, m and l are the number of j 

4' 

parameters represented by the row vectors E, F and G respectively. Simi- 5 

larly, the result of equation (2. 59) is a 3 s 1 matrix. These equations are j 

to be used together with equation (2. 25) to find the partials of the measured 

I 

distance d t with respect to these physical librstlon parameters. ; 

The GAST I 

■■ r 

The Greenwich Apparent Sidereal Time (6) is given by (see [69]): 

j 

6=-- UT+ 6 h 38 B 45S836 + 8, 640, 184/542 t„ + 0!0929tf + A# cos e (2.60) j 


where t^ is the number of Julian centuries of 36525 mean solar days since 
1900 Jan. 0. 5UT. UT (Universal Time) is the UT1 epoch of observation, 
the next three terms gives the Mean Sidereal Time (MST) at C^UT, and 
the last term is generally referred to as the equation of the equinox. A$ is 
the nutation in longitude and e is the true obliquity of the equator of date. 

Each UT1 time of observatic •• cannot be taken as a parameter, since 
this will always introduce one more parameter to the system each time an 
observation is made. However, general time corrections (clock offset error) 
cab be found for fairly short spans of time. If UT1* is the chronometer 
time e£ observation, and AT is the "regional” correction to UT1* in order to 
obtain the correct UT1, then 

UT1 = UTl* + AT. 


Together with equation (2. 60) 






If equation (2. 60) is rewritten as 


where 


© = UT1 + I?T + A(J> cos €, 


= T t 


T = t„ , and B = B ; 


It follows therefore, that 


ax T _ ax, ae. _ ax, 
a&r a© * a&T a© 


(2. 61) 


ax T _ ax, a© 
aB a© ’ aB 


(2.62) 


iu given by equation (2. 41), and ■ can be obtained using equation 

a© diiT.oiJ 

(2.25) as demonstrated earlier. 

Since the "variable" parameters and c are also used in computing 
the value of ©, equations (2. 51) and (2. 46) can be modified by adding 
respectively the expressions: 


ax T ax, # _a©_ d£*b 

ax a© ’ a&£ ‘ a© 

ax T = . . + ax, . a© # a? 
ar%/'»a3] T a© at arao»’**<ibl 


(2. 63) 


(2. 64) 




The first part of equations (2.63) and (2. 64) are given by equations (2.51) 
and (2. 46). From equation (2. 60) 


= cos ( 


= -Atfi sin c. 


In practice, since cos e is a small quantity, the equation of the equinox 
contribution to 6 can be regarded as errorless and if this is so, it i3 not 
necessary to use equations (2.63) and (2.64) in place of equations (2.51) 
and (2. 46) respectiv >ly. 

Coordinates of the True Pole 

The coordinates x P , y P of the true celestial pole (with reference to the CIO 
pole) are traditionally determined from the analysis of continuous latitude 
and/or longitude observations at permanent observatories. Two agencies — 
the International Polar Motion Service (TPMS) and the Bureau Internationa! 
de l'Heure (BIH)-use the variation of latitude and/or longitude values at these 
observatories in determining the motion of the pole. The IPMS uses primarily 
five stations located near a single parallel of latitude, while the other stations 
participate primarily in the BIH program. 

The coordinates of the true pole are published by both the IPMS and the 
BIH. The final coordinates are given by IPMS at intervals of 0. 05 year 
(18. 25 days) but published two years in arreas. The BIH final coordinates 
are given at intervals of five days and are generally available one month 
in arrears. 

To date, the theory of polar motion is only approximately known. 

There are no existing mathematical functions which can be used to calculate 
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the position of the true pole at any epoch such as those that have already 
been given for other parameters. However, coordinates of the true pole — 
x p , y p — cannot be regarded directly as parameters as was noted above. 
Hence, corrections to known values of Xp, y P will have to be assumed con- 
stant for a short length of time (such as five days), or alternatively, the 
corrections can be obtained in another way. 

Recently, polar motion have been obtained from residuals of Doppler 
satellite observations by J. Anderle and Beuglass of the U.S. Naval Weapons 
Laboratory [ 8 ]. In a similar manner, corrections to the values of Xp,y p 
used in predicting laser distances can be obtained by analysing the residuals 
of the measured distances after adjustments. 

The Geocentric Position of the Selenocenter 

The position of the moon's center of mass given by Xcq, Ycq, Zc q are, 
like other variable parameters considered above, computed for each epoch. 
For this purpose, a lunar ephemeris, analytical or numerical, has to be 
used. The lunar ephemeris with the smallest estimated uncertainties is 
the JPL LE16 which is estimated to have uncertainties of 100-150 meters 
in position and 50 meters in range. The laser ranges are expected to be 
at least one or two orders of magnitude better in precision than the best 
ephemeris. It is therefore improper to use any of the existing lunar 
ephemerides in the prediction equation without allowing for the systematic 
errors in the ephemeris used. 

In order to take care of the systematic errors in the existing lunar 
ephemeris, such as the numerically integrat . LE16, an empirical model 
can be used in the observation equations, who* - sole purpose would be to 
absorb the systematic errors in the ephemeris. This will not necessarily 
lead to an "improved” lunar ephemeris, but will prevent the distortion of 
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the adjustment due to systematic errors of the ephemeris. Such an 
empirical model can be a harmonic series, with fixed coefficients and 
arguments which are linear functions of time, or a general polynominal, 
the degree of which will depend on some pre-set conditions on the analysis 
of variance obtained in the adjustment. 

Such a harmonic series can be of the form: 


£ fit cos (3|t + b,). 

i 

Thus, the "true" geocentric coordinates of the selenocenter that will be 
used in the distance prediction equation as well as in the adjustment model 


can be of the form: 


XqT] £ fj^ cos(0,‘t + b,*> 

t 

Yqc + Z fl? COS (fit t + b, v ) . 
i 

Zg. £ ft} co3(fi} t + b,*) 


If a polynomial is used to model the systematic errors in the ephemeris, 
the "true" geocentric coordinates of the moon will become 

fXqcl fco + c x ‘t + Ca t 3 + .^1 


Co + C x Y t + C 3 Y t a + ••• 
C^ + C't + Ca t 3 + 


Depending on which model is used, the coefficients p,, /3 t , b, or C t become 
part of the parameters to be solved for in the adjustment. 

The Geodetic Coordinates 


Equation (2.44) gives the partiais of the topocentric coordinates of the 
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moon point with respect to the U, V,W (geodetic) Cartesian coordinates of 
the laser station. The U,V,W coordinates are functions of the geodetic 
coordinates (o, X, h) of the earth point and the parameters of '.he ellipsoid 
(a, f) (see equation (2. 1)). Optionally, if the position of the earth point is 
defined in an absolute geodetic system, one can regard o, X, h, a, f as param- 
eters in place of the Cartesian coordinates. Consequently, equation (2.44) will be 
replaced by the following equations: 


ax, _ 


(M + h) sino cos X 


= -R^CeJPVs 1 -(M + h) sintp sinX 
(M + h) cos <p 


(2. 65) 


-(N + h) coso slnX 
It 1 = -Ri(€o)P t N t S t (N + h) cosocosX 

O A 


( 2 . 66 ) 


cos o cos X 


= -R l (c 0 )P T N T S T cos O sin X 


(2. 67) 


= -R^foJPVs 


cos o cos X 
W 

iTj/gT COSO sin X 

(l-e 8 ) sin o 
W 


( 2 . 68 ) 




M = 

W = (l-e a sin^pf^ 

and a, e, f are raspestively the semimajor axis, eccentricity and flattening 
of the ellipsoid. 

The above equations (2. 65) to (2. 69) are valid only when the U, V,W 
coordinates are expressed in the "absolute" geodetic ("average" terrestrial) 
coordinate system. In many instances, coordinates of earth stations are 
expressed in a "relative" geodetic coordinate system, using any of the 
conventional reference ellipsoids as the reference figure for the earth. The 
reference ellipsoid used may not be centered at the geocenter, and the 
Cartesian axes of such a "relative" system may not be oriented parallel to 
the axes of the "average" terrestrial coordinate system. If the relative geodetic 
(curviilinear) coordinates of the earth station is ip, TL, and E, and the param- 
eters of the reference ellipsoid are a and then the "relative" Cartesian 
coordinates u, v,w are also given by equation (2. 1) provided <p, X, h, N and e 
are substituted for <p, X, h N and e respectively. A general relationship 
between the "average" terrestrial and the "relative" geodetic coordinates of 
a point are given by (see [68]): 


afl-e 8 ) 

W 3 
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V 1 


U u] 6 0o 6 A 0 A u 

V - vj+A6X 0 + Bfi£o+s&v 


( 2 . 70 ) 


U - Uq 


A v = v - v. 


! A w 


sin Oo cos Xo sin To -cos 0 O cos To 
sin 0 O sin To -cos To -cos0 o sin 
jcos 0o 0 -sin 0o _ 


! -cos0qAv- sin0oSinXcAw 
, cos0bAu + sin0bcosXoAw 


cos0bsinXoAu-cos<£bcosAoAv cosXoAu + sinXoAv sintOosinAoAu- simcbcosXoAv 


sin<Cb Av - cos03 sin Xo Aw J cos <Oo Aw 

i 

B = -sin^j Au + cos0jcos X^j Aw < sinXoAw 


00 — 00 


6 X<3 = Xo — Xo 

5 ho ho - ho 


6 Ao is the error in azimuth of the relative system 
6 Co Is the error in tilt of the system in the meridian plane 
6 rjo is the error in tilt in the prime vertical plane 


s is a scale factor. 


It should also be noted that 60o and 6Xo should be expressed in linear units. 
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All the variables with zero subscripts denote the values of these variables 
fit the origin of the relative geodetic system. 

If equation (2.70) is now substituted in any of the equations (2.10), 
the topocentric coordinates of the moon point (X T ) will be obtained. Also, 
in equation (2.70), seven transformation parameters (&Ao, 6X0, Sho, 6Ao, 6£o» 
6r)o and s) have been introduced in addition to the "relative" geodetic coordi- 
nates and the parameters of the reference ellipsoid. The partials of X T with 
respect to ip, X, E, a andT are given by equations (2. 65) to (2. 69), provided 
the relative quantities (with bar) are substituted for their absolute counter- 
parts. Taking the partials of X T with respect to the seven transformation 
parameters we obtain: 

= -Ri(Co)P T N T S T A (2.71) 



and 



-R l (fo)P T N T S T B 




(2. 72) 


(2.73) 


The above three equations cau be used with equation (2. 65) to (2. 69) 
instead of equation (2.44) if it is desired to obtain corrections to the relative 
geodetic coordinates of the earth station, the parameters of the (relative) 
reference ellipaoia as well as the seven transformation parameters. 
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The Selenodetic Coordinates 

-■.s in the case of the geodetic coordinates, the Cartesian coordinates of 
the moon point (x«, y», z*) ari functions of the spherical coordinates (t, b, r) 
of the lunar point (see equation (2. 2)). Equations for obtaining partlals 
for the x*. y M , z M coordinates have been given in equation (2. 43). If, however, 
the spherical coordinates b, l, r replace x„, y M and z„ as parameters, then 
the following equations replace equation (2. 43): 


= Ri(fo)P T RiH)Tl 


-r sinb cost 
-r sinb sint 
r cosb 


(2.74) 


- Rx«o)P T RiH)Ti T 


-r cosb sin-t 
r cosb cost j 
0 


(2. 75) 


ax r 

dr 




cosb cost 
cosb sint 
sinb 


(2.76) 


I 2. 5 Observation Equations when Orientation Angles 

i are Obtained t * rough Numerical Integration Process 

| The idea of obtaining the moon's Eulertan angles through numerical 

i 

■ integration of the equations of motion have been discussed in Section 2. 34. 

1 Also, it has been suggested in Section 2. 35 that the conventional precession. 

I nutation and the Greenwich hour angle of the vernal equinox can be replaced 

by three (Eulerian) angles, which can also be obtained by numerical inte- 
\ gration. The numerical integration procedure for obtaining these earth 

i 
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and 


E* - • ip* ; 4"*! 

respectively represent the orientation angles of the earth's UVW and the 
moon's xyz axes (fixed to the bodies) with respect to an inertial system (such 
as die mean ecliptic coordinate system of 1950. 0), then the topocentric 



coordinates of the moon point is given by equation (2. 15): 


x* 


| ^4; 

h 

c" 1 

i 

Y„ 

= Ri(co) 

! Y cp ! 

+ p« i y» 

i 

- P E V 

z„ 

U J 

topoe. 

: Zcq 

W mJ 

L Zm 

w 


where P t are 3*3 orthogonal matrices given by: 

P« = ^W>*)Ri(e«)Ifcf-&K) 

P E » lM-<MRa<9 E )R*H’£) - 

The predicted distance from an earth station to a moon point is computed from 



equation (2. 16): 

D = P& + Yg, + Z^)*. 

The matrix of partials (design matrix A) will be different from those 
already derived in Section 2. 4, since the parameters are different from 
those listed in equation (2.18). In analogy to equation (2. 18) 

= f3(£ot$’H>Qn»4‘n»!{ , E »9 E « V,W, x^y^, z Ht Xcq> Ycq, Zcq) (2.78) 
In addition, since the Eulerian angles are obtained from integration. 
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The variable parameters (0„, • • • , <J> E ) will have to be replaced by the 12 
initial conditions in equations (2. 79) and (2. 80). If 




ESI = 0?2 0 2 *2 fi 2 <*2 $21 


Then 


dX T _ BXj BE„ 
3[eS : ES] be*' cfES ; e2i 


(2.90) 


where — J is a 3 x 3 matrix of partials whose columns are given by equations 
dE M 

(2. 82) to (2. 84). Similarly for [E° : Ef», 


BX, _ BX T BEr 
BfE° l Ep BE t * 3[E° i t%] ' 


(2-91) 


Also, — 1 is a 3 * 3 matrix of partials, whose columns are given by equations 

oEg 

(2. 86) to (2. 88). Each of the expressions ...a and -.r^O ■ * ‘fr. i8 a 

: EbJ : t-cJ 

3*6 matrix which forms part of a matrix, usually referred to as the state 
transition matrix T 9 ]. This matrix describes the transition of a differential 
variation of the initial epoch conditions from the initial epoch t 0 , to time t. 

The matrix can be obtained along with the integrated orientation angles by 
methods which are described in Chapter 4. 

The partials of the laser distances with respect to each of the parameters 
considered in this section can be obtained in the same way as was done in 
Section 2. 4 through die relationship given by equation (2. 23) 

BD _ BD BXt 
bx bx t a» 


where D represent the distances, X T moon point tcpocentric coordinates and 
x the parameters. 
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2. 6 Summary 

In this chapter, one of the new observational systems which coaid be used 
to improve selenodetic control has been considered. A brief description of 
the laser ranging systems was given and equations relating ranges between 
an earth point and a moon point to astronomical, geodetic and selenodetic 
parameters were derived. The range equations are valid for laser ranges 
as well as any other kind of ranges made between a station on the earth and 
a moon point. 

The range equations were differentiated for the formation of observation 
equations which can be used with lunar range measurements in an adjustment 
model to obtain corrections to approximate values of the parameters involved. 

In deriving the observation equations, observed ranges are assumed to 
be corrected for atmospheric, instrumental and other systematic errors. 

In Appendix B, the effect of the atmosphere on observed laser ranges is 
briefly discussed. 

The range equations can be used in predicting laser lunar ranges for 
real observations, and can also be used to simulate observational data for 
experimental adjustments as was done in Chapter 5 of this work. 
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3 APPLICATION OF THE VL3I TO SELENODETIC CONTROL 
3. 1 Introduction. 

Radio Interferometry began in 1346, when McCready et al., used an inter- 
ferometer for solar observation [17], Later in the same year, the 
two-element interferometer also came into use at Cambridge. Radio 
interferometry rapidly developed thereafter with longer baseline obser- 
vations, but the maximum length of baseline was limited. The conventional 
interferometry used superheterodyne receivers at each terminal. A 
common local oscillator signal had to be transported across the baseline 
by means of a cable connecting the two terminals. Over longer distances, 
the operations were managed across the baseline by using a radio link 
containing two repeater stations. The length of baselines for conventional 
Interferometry hardly exceeded 130 km. 

The development of atomic frequency and time standards found a 
useful application In radio interferometry. Using independeut atomic 
standard oscillators, all real-time interconnections between the two inter- 
ferometer stations could be eliminated. The interferometer signals from 
each telescope are independently recorded on a magnetic tape along with 
time, and the tapes are later brought together for processing. This tech- 
nique of "atomic-clock" interferometry or "independent clock" interferometry 
(with local freqency standards and tape recorders replacing microwave links 
or cable intercoauections) is now known as Very Long Baseline Interferometry 
(VLBI). 
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With the VLBI, resolution of distant radio sources with a n gula r 
diameter of0"0006 has been achieved. This compares with a resolution 
of ^ minute of arc achievable by the largest reflectors, which is roughly the 
same 3s the limit for the unaided human eye [16], The angular resolving power 
of a radio interferometer system is given by 



(3.1) 


where 


a is the angular resolution in radians 


and 


A is the wavelength at which the antennas operate 
d is the length of the baseline. 


It is therefore necessary to know the length of the interferometer baseline 
in order to determine the size or position of the radio sources. 

The inverse problem is of more importance to geodesists. If the 
coordinates of point radio sources are available, then the baseline param- 
eters of the interferometer system can be determined. 

Many possible applications of the VLBI in geodesy and geophysics have 
been suggested by various scientists in recent years. Such areas of appli- 
cation include geodetic ties between continents, continental drift, deter- 
mination of precession, nutation and the earth's rotation rate, satellite 
tracking, navigation and time synchronization. In this chapter, the 
possibility of using the VLBI for determining station positions on the earth 
and on the moon will be investigated. The relative position of a lunar 
station to an earth station at any instant is a mathematical function of other 
parameters apart from the geocentric and selenocentric positions of the 
earth station and the lunar station respectively. Consequently an updated 
knowledge of these parameters can be expected when the VLBI is used 
for earth and moon position determination. 
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It is envisaged that the establishment of radio antennas on the moon 
in the future is a possibility, notwithstanding present technological limita- 
tions, Consequently, the studies in this chapter Include Interferometric 
systems made up of earth-to-earth, moon-to-moon and earth-to-moon 
baselines. 

In the next section, the basic operating principles of an interferometer 
system will be given. An attempt has been made to make this brief descrip- 
tion as non-technical as possible. Its presentation in this chapter should 
point out the type of observations to be expected from the VLBI, thereby 
making the rest of this chapter more easily understood. 

3.2 Basic Operating Principles of the VLBI. 


An Interferometer is an apparatus for measuring the phase difference 
between simultaneously received electromagnetic radiation at two stations 
from a distant source. Alternatively, the elapsed time between 
the arrival of any particular wave front at each of the two stations (time 
delay) can be measured. The phase difference & p and the time delay r 
are related by the equation 

r = 4*- (3.2) 


where 


u) is the nominal frequency of the radiation. 


The necessary equipment at each of the interferometer stations include 
a radio antenna dish, a local oscillator and an atomic clock (both controlled 
by a highly stable frequency standard such as a hydrogen maser frequency 
standard generator). Also there should be at each station a mixer, a video 
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converter, clipper, sampler and a tape recorder. These pieces of equip- 
ment are schematically shown in Figure 3. 1. 


Terminal A 


Terminal B 



Figure 3.1 

Equipment for a VLBI Observing System 
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The look-angle of the radio source (such as the quasar) at the desired 
time of observation is pre-computed at each station and each radio antenna 
is steered to point toward the direction of the radio source. The receiver 
is also tuned to the desired operating frequency. A synchronizing circuit 
starts the recording of radio signals received and time, at a pre-set 
epoch. 

The recetved signal, whose central frequency is w, is mixed with a 
local-oscillator produced frequency that is phase-locked. The resulting 
frequency out of the mixer, known as the intermediate frequency (IF) is 
passed on to a video converter. The video converter’s job is to convert the 
output signal into a low-frequency signal, suitable for 1-bit digital sampling 
and recording. The IF signal is amplified, filtered and phase-coherently 
converted in frequency to obtain a single sideband output signal whose 
carrier frequency is zero. This output video band signal (limited to 480kHz by 
available digital recorders) is then passed on to be clipped, sampled and 
recorded. 

Clipping. 

The received signal now converted to video band frequency is com- 
pletely random, a characteristic of the emitting natural radio sources such as the 
quasars. This signal is clipped "infinitely ", thereby retaining the sign of 
the voltage received and dispensing with the magnitude (see Figure 3.2). 

A correction can ’a ter be applied to the correlator output to compensate for 
the error introduced by clipping the original wave form. The correction 
to the correlation output according to the Van Vleck clipping theorem 
is computed as, [84]: 

Mt) = si" [f R,(t>] (3.3) 
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where 

Rx(t) is the corrected signal 

R T (t) is the raw output of the correlator 



— > time — > t 


received signal clipped signal 

(a) 

Figure 3. 2 

Clipped Signal in Relation to the Received Signal 

Sampling . 

For the purpose of recording the clipped signal on a magnetic tape, 
the sign of the signal is identified at a regular, spaced interval of time. 

The sampling rate 's, in accordance with the Nyquist theorem, equal to 
twice the recording bandwith. Thus as an example, the National Radio 
Astronomy Observatory (NRAO), which uses a bandwidth of 360 kHz, records 
the data at 720 Kblts/sec. In addition to the clipped and sampled signal. 
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time is also recorded at intervals of 0.2 sec. The recording of signals 
and time can be done either on a video tape or on a digital tape. It is 
more convenient to record on the video tape when sources are be ug ob- 
served since a video tape can last many hours of observations, whereas 
a digital tape can only record three to four minutes of data. If video tape 
is used, it is necessary to re-record the data on a digital tape afterwards, 
for digital computer use. 

Data Processing. 

The data now on the digital tape consists of the clipped and sampled 
signal with time, one tape for each of the two interferometer stations. The 
processing of the tapes is done in two phases, namely: 


(a) correlating the tapes 

(b) obtaining the true delay. 

The correlation of the tapes is done by a digital computer. For each record, 
the bit streams are shifted relative to one another by an amount that corre- 
sponds to the theoretically predicted geometric time delay. The predicted 
time delay is a function of the baseline length d , and the angle 8 which the 
baseline makes with the direction of the observed radio source, 

T « 

where c is the velocity of light, 
given by 


d cos g 
c 


(3.4) 


The theoretical correlation function is 


R(T) 


llm 

T-» 00 


1 _ 

2T 


£ x(t)- x(t + T)dt 


(3.5) 
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Figure 3.3 Basic Geometrical Relationships in Interferometry. 
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Clipped signal received at A 
as function of time, starting 

aft,- 



Clipped signal received at B, 
starling at l- . Note the shift of 
signal at B v. r. to A 


Figure 3.4 Clipped Signal Received at Two Stations and Time Delay, 
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In practice because of finite record length the function is given by numerically 
integrating over N records, i.e. 

"1 

R'(t) = £ x(to + nAt) y(to+n6t + r) (3.6) 

B— l 

using a number of values for r centered around rhe predicted geometric 
delay r g , and N t (number of records over which integration is to be carried 
out). N t depends on many factors such as 

(1) source strength 

(2) antenna size 

(3) system noise of receiver 

(4) clock stability, i.e. phase stability of the frequency standard. 

The first three items put a 'ower limit on the desirable integration time since 
weaker sources demand more integration time, and so also do small antennas, 
and receivers with high system noise. On the other hand, stability of the 
clock used places a higher limit on the integration time because phase stability 
must be maintained to high accuracy (1 part in 10 a or better) within the inte- 
gration interval. 

The correlation function obtained is corrected for finite record length by 
multiplying R'(t) by a weighting function w (r) which is chosen to satisfy the 
following criteria: 

w(0) = 1 (normality) 
w(-T) = w(r) (symmetry) 

W < T >T2T B = 0 

where T n is the total length of the record. The last condition ensures that the 
corrected function r'(t) is defined for all r. An example of such a weighting 
function used in practice is given by the ,r Hanning" function r 843: 
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W(T) = J[1 + COS ~J, if |T| <T 0 


W(T) = 


, if |r|>T.. 


The weighted correlation function is plotted against time delays and a 
sharp peak in the graph is an indication of good fringes. 

In the plot of the correlation function versus time delays, the central 
peak of the plot denotes the approximate position of the true delay on the 
abscissa (see Figure 3.5). The exact peak is determined by fitting a 
theoretical function to the plot of the correlation function. In a perfect 
situation, the observed true delay (r,) should be the same as the predicted 
geometric delay (r s ). Any difference between these two time delays (t b - ij) 
is due to errors in apriori knowledge of baseline parameters and other effects 
such as clock offset error, differential refraction and instrumental delays. 


correlation function plot 

cu i vc fit 

r t - calculated delay 
7, - observed delay. 


! ;h 


r r t T » t,i4t 

Figure 3.3 Fringe Amplitude Plotted Atrainst Time Delays. 

The accuracy with which the time delay is measured depends on the receiver 
bandwidth, among other factors. When a receiver is teaed to a central frequency <4 
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it receives a broad range of frequencies centered about w. The width of 
this range depends on the receiver characteristics and if this bandwidth is 
Aw, then the accuracy of the delay measurement depends on the correlation 
interval At given by 


At 


1 

2- Aw ' 


(3.7) 


Tliur the accuracy can be Increased by narrowing the correlation Interval, 
i. e. Increasing the bandwidth. Although there are practical limitations on 
the bandwidth attainable with a receiver, the effective bandwidth can be 
made larger by sampling narrow band video signals from many widely 

separated windows. This method is described in detail by Hinteregger [34]. 



It has been shown in the last section how the VLBI operates and that the 
observed quantities of interest to geodesists is the difference in the time of 
arrival at two separated antennas of a particular wavefront from a given 
point source of radio radiation. This relative time delay (or phase delay) 
is a function of the source direction, the antenna locations, the relative 
clock error between the two sites and any other influences such as the 
differential atmospheric refraction. By expressing the measured time delay 
as a function of the various parameters involved, and with sufficiently large 
number of observations for each of the observed celestial radio sources, the 
"best fit" values of the parameters can be estimated by weighted least- 
squares analysis. 

In this section, expressions relating VLBI measured time delays to 
parameters which Influence the measurements will be derived for different 







configurations of interferometer stations. The three basic configurations 
envisaged are those VLBI observations involvbg: 

(1) two stations on tne earth (earth-earth VLBI) 

(2) two stations on the moon (moon-moon VLBI) 

(3) one station each on the earth and on the moon (earth-moon VLBI). 

Although the primary interest in this study is the investigation of applica- 
tion of VLBI to selenodetic control, other parameters indirectly connected with 
the determination of coordinates of points on the moon will also be considered. 
Thus, the case of the earth-earth VLBI will be treated even though selenodetic 
control cannot be derived directly from the observations except in the case of 
the earth stations observing an artificial radio source on the moon. 











Two radio antennas A and B situated on the earth have position vectors 
-> — ► 

p A and ft respectively. The inter-site vector D is given by 

D = ft - p A . (3.8) 

-4 

The position vectors of o* and ft can be obtained from the geodetic coordinates 
of A h*) and B&fe, X a . M as 

p A = U A r + VaT + W A k 

^ (3. 9) 

ft = Ub? + V B T + W a k 


where 


"u»l f(N,+ 

Va = i (N a + 


1“(N» + h*) cos<& cos X A 
| (N a + h*) cos <& sin X A 


a J |_(N» (l-e 3 ) + h*) sin 


(3.10) 


(1^ + he) cos (Da cos X B 
(hfe + he) cos cfe sin X s 


t — of — 

_(*k (l-e 3 ) + 


he) sin (ft 


(3.11) 


D = 4ji + d»J+ d^k 


(3. 12) 


where 


(l-e a sin 2 o)® * 

^ 4 4 

Also, i, j,k are unit vectors along the U,V,W earth-fixed coordinate system. 
Thus the vector D expressed in the U, V,W coordinate system is given by 




."Si s. 


( 3 . 13 ) 




Let 3 represent the unit vector from the radio source to the geocenter 
Natural radio sources such as quasars are nt extragalactic distances away, 
so that the geocentric parallax can be neglected. Thus, s also represents 
the vector from the source to any point on the earth. The vector sis 
obtainable from the catalog of radio sources which gives the celestial posi- 
tions of the sources in a celestial coordinate system (usually the 1950. 0 
mean equatorial system). If the cataloged source position is ckq, 6 0 » then the 
components of the unit vector Sq expressed in the 1950. 0 mean equatorial 
Cartesian coordinate system is given by the direction numbers as 


So 


-cos cos ol 
- cos 6 o sin 

L "Sin 


( 3 . 14 ) 


The true equatorial Cartesian coordinates of the source at the observation 
epoch is obtained from So by applying precession and nutation. 


- i. 
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The delay T is a function of the baseline length d and the source direction 


with respect to the baseiine as 


d cos 


, c ~ light velocity. 


Thus, if the vectors D and s arc expressed in the same system, then since 


s is a unit vector. 


1 -» -* 

T = ~(D. s). 


(3.15) 


The above equation would give the difference between the two times of 
arrival (tg. t : ) of the wavefront at B and A if station B remained stationary 
during the interval (t^— t t ). However, due to the change in the position of 
B in the interval (ts-t^, the true delay is 


' 5-1 


(3.16) 


where pa is the velocity of site B in the same coordinate system in which 

— > -p 

s and D are expressed. 

The vectors D and s are two vectors in space, and their scalar (dot) 
product would remain the same irrespective of the coordinate system in 
which the components of the two vectors arc expressed. It is only impor- 
tant that the components of vectors be expressed in the same coordinate 
system. For convenience, as well as for the desire to work in an inertial 
coordinate system, the 1950.0 mean equatorial coordinate system can be 
chosen. Thus, if s 0 and D Q arc the vectors D and s whose components are 
expressed in the 1950.0 mean equatorial coordinate system, then 


Dr, * rT] D 


(3. 17) 


SO 





and Sq is given by equation (3. 14). [T] is a transformation matrix necessary 
-» •* 
to transform the D vector expressed in the UVW coordinate system to D 0 

vector expressed in the 1950.0 mean equatorial system (XYZ). Thus 

m = p i n i s i 

where P, N, S are 3*3 orthogonal transformation matrices, whose full express- 
ions have been given earlier in Chapter 2. Equation (3. 17) can now be written as 




Do = [P t N t S t ]D. (3.18) 

In a similar manner, the vector ft (in equation (3. 9) can be expressed 
in the 1950.0 mean equatorial coordinate system: 

ft 0 = fP T N T S T 3ft. (3.19) 

In equation (3. 16), the rate of change of the vector Pb 0 with respect to 
time is required. Rewriting equation (3. 19) by substituting expressions 
for S T , N t and P\ 

+ AORa(-®)Ri(yp)Ra(x P )]p 6 . (3.20) 
Differentiating equation (3. 20) with respect to time: 

* r<&,LoP T N T S T ) -<fy*a(W La R 3 Hi)R3<Zi)N t S t ) + 

+ (z l P T L 0 N T S T ) - £ P T Li N t S' ) + (€P t N t LiS t ) + 

+ (6$ P T R l ^)L Q R a (^)R^ + Ac)S T ) + (A<P t N , L,S i ) + 

- (6p t N t LoS t ) + (y ( ,P I N T R 3 (-e)L l R l (yp)R s (x P ) + 


+ (ic, P T N 7 S T L^1 ft. (3.21) 

In the expression for ft Q , the subscripted L matrices are the Lucas matrices 
defined in Chapter 2. The dotted quantities are the rate of change of those 
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quantities which can be obtained by differentiating the expression for each 
quantity with respect to time. Analytical expressions as functions of time 
do exist for each of the dotted quantities except for X? and y p , whose values, 
though . small, can be calculated from tabulated x* and y f values. 

The pa vector, expressed in a geodetic UVW coordinate system, is invariant 
with time. 

If the vectors Dt» ’k> and £a 0 have been evaluated from their equations 
given above, let 

Db = Djf + Dj,? + Dak 

sb = s, ? + s a ? + 83 k (3.22) 

k Q = V t i + v a j + v 3 k 


where 

-> *♦ -ft 

i, ], k now represent unit vectors along the X,Y,Z axes that 
define the 1950. 0 mean equatorial Cartesian coordinate 
system. 


Then equation (3. 16) can be written as 


c-r = (DiS! + Efesa+ Ebsi,) 


Equation (3. 23) represents the prediction equation for the measured 
time delay for VLBI observations involving two stations which are located 
on the earth's surface. The predicted values are compared with the 
observed values and the differences can be used in an adjustment process 
to correct the parameters which implicitly appear on the right hand side of 
equation <3. 23). 

Alternate expressions for the equations for the two vectors D 0 and 
3 bo can be derived by considering that the orientation of the UVW coordinate 
system with respect to an inertial system is equivalent to three Eulerian 




j 

I 






(orientation) angles at observation epoch. This concept, and how to obtain 

y 

the orientation angles and their time rates i-= given in Chapter 4 of this work. yj 

Let 9, 0 and <it define the orientation of the UVW coordinate system ,| 

with respect to the mean ecliptic coordinate svstem of 1950.0. Then, the 

vector D 0 is given bv f. 


D 0 = r «.(-Ca)T T lD = rK.f-C'JK.H)K : (^K>H4lD (5.24) 

where 

T - Ui(4)U,(-9,K < c) 

Co is the mean obliquity of the equator at 1950.0. 

Similarly, 

Ofi 0 ^ rR 1 (-Co)H3e^«.(9>Ha(-4)iP3 = r rti(Co)T r lft (3.25) 

Differentiating equation (3,25) with respect to time 

= [-(0Ri(-Co)L a T T )+( < 5R 1 (fo)R3K)L l H l (e)R3e-H)U-^Ri( f o)T T L^]^. (3.26) 

Since ty, 9 and 4 « are obtained from numerical integration as well as 9, 0 and 
the equations for Do vector and vector can be evaluated. The predicted 
time delay can then be evaluated from equation (;;.23). 


Observation Equations for Earth- Earth VLDI . 

In order to use the differences between the predicted and the observed 
time delays in improving the present knowledge of the parameters involved in 
VLB! observations, it is necessary to derive the partial derivatives of the 
time delay with respect to the parameters. The paramci s involved arc 
the station coordinates, the positions of the sources, t jjreeessional ele- 
ments, nutation and polar motion parameters and the -'"vieh apparent 
sidereal time. Alternatively, the Eulerian angles 9. C an be regn>- 
as parameters, in addition to the positions of the statim. and rad i sources. 

The equation for the predicted time delay (equation (3,23)) can be written. 


i 



! \ 


v : 4 

l " I 

' rf 


! . 


In matrix notation, as 

T = ~ T(D t s) + A,(D t s v T 8) 

where 


D 0 = D 


D, 

Ds 

D, 


So = S 


% 


and 


Vl 

V = I Vg 
V 3 _ 

The parttals of r in equation (3. 27) are given hy 


TT = ~ r S T l + fs T V s r ] 
oD C C 


I “ - ~CD T ] + pftDW) + (v t 8D t )1 


dr _ _l 

dv e 


= ~ a F s T D a T ] 


84 


(3.27) 


(3.28) 


(3. 29) 


(3. 30) 


.'4 


r A 


I ■ '.■] 

\ I 


■m 


Vi f! 
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Each of the above eouations for partials results in a 1 * 3 row matrix. The 
elements of the matrices D, s, v are functions of parameters, and the 
matrices will have to be differentiated with respect to these parameters. 
Hence if the parameter list is represented by vector X, i.e. 

D f,(x) 


v = (x) . 


bT _ OT nj) | bT r)S ^ OT bV 

bX -)D bX bS bX bV bX 


(3.31) 


where the partials — , — and ~ arc obtained bv differentiating expressions 
bX OX ox 

for D. s, and v. 


From equation (3. 14) 


-cos fy 0 cos a 0 I 


s = -cos 5 0 sin a 0 . 
-sin 6 0 


Hence 


sin fj 0 cos a 0 

bs ... 

— = sin 5o sin a 0 
OO 0 

-COS 5a 


(3. 32) 


(3. 33) 


ds 

too 


cos 6 0 sin »o 
-cos h 0 cos at# 
0 


so that 4t . 

O0q 


dr 

da 0 


are obtained using equation (3.31). 


From equations (3.18) and (3.13), 


D = (IV>.(Cu)Ra(-8i)Ka(3'.,)R l (-OK n (^|i<i<c + Ac)*W-e>Ri(y»)Ra (*>)) 


Hence 


dD 

oCo 


L 0 P , N t S t 


d. 

d. 


d* 


“ -lh(r o) La R a ee,)R,(it 1 )N T S T 



dD 

07.! 


P’LjN’S 1 


dj 

dy 

d* 




~ d. 


"d. * 

oD 

- -p t l,n t s t 

d, 

+ P'.VX.S’ 

d, 


_ d «_ 


_dy,_ 


cL 

dy . 
d w 


(3. 34) 


(3. 35) 


(3.36) 


(3. 37) 
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dD 

d&ii> 


P'R^OIqI^AMRi* +AC)S' 


<3. 38) 


4j 

d, 

d« 


dD 

d&e 


P t N t LiS t 


d, 

d; 


d* 


dP 

d© 


= -pVloSf 


d. 

dv 

d* 


(3. 39) 


(3.40) 


dP 

dy» 




4 

dv 

d* 


(3.41) 


Since 


dP 

dx P 




P T N T S T La 



X ■ 


Ub" 

ku 

dv 

= 

v 8 

-!v. 

> 


_W e _ 

k 


(3. 42) 
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-P r N r S r (3.43a) 


P T N t S t • (3.43b) 


As can be expected, both o» and pg vectors cannot be solved for together and 
either of the vectors have to be fixed. Also, all the parameters whose 
partials are given by equations (3.35) to (3.42) are ''variable'' parameters, 
and their values vary from instant to instant. However, each of the parameters 
can be expressed as functions of another set of "non-varying" parameters. 

It is desirable to replace the 'Variable" parameters with the "non- varying" set, 
and this can be done with the same equations given in Section 2. 4 of this 
work. 

The alternate expression for D, from equation (3.24) is given by: 



'<*»“ 


‘clT 


d, 

- d "_ 

= R:(-eo)T T 

d, 

d* 


so that 


3D 


U* 

V* 


3D 


3D 

3c 0 


-L, R 1 (-c c )T r 



(3.44) 
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7? "l 
* A 


bP 

dip 


= -RtKo)I« T T 


4, 

d, 

d* 


(3.45) 


m- 

*g 

'I 

i 

■ % 


ap 

30 


R ! e€o)i\sH)LiR 1 (6)R 3 (- < i) 


4, 

dv 

dn 


(3.46) 


| % = -R^€o)T t Lo 


4j 

dv 

d» 


(3.47) 


3D 

U» 

V* 


W u 


= -R,K<JT T 


(3.48) 


and 


= R,Ko)T\ 


(3.49) 


3D 
Ua 

V 8 
L W 8J 


The "variable" parameters 0, 0 and 4- can be replaced by the initial conditions 
(0o. ii)o. * 0 , 6o, Vo and <i>o) and the partials of r with respect to these "non- 
varying" parameters can be obtained by using the state transition matrix 
(see Chapter 4), equations (3.45) to (3.47) and equation (3.31). 
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The expression for p& 0 given by equation (3.21) should also be differ- 
entiated with respect to the parameters in order to define all the terms 
appearing in equation (3. 31) (the v matrix is represented by the vector pe^. 
The largest of the terms in equation (3.21) is that dependent on the sidereal 
rate, t. e. , 



v ^ -©(PVLaS 7 ) 


U b 

V a . 
W 8 _ 


(3. 50) 


Partials of equation (3.50) are given by: 



-(PVl«S T ) 


U b 

V b 

_W*_ 


dv 
a Co 


-©(LaPVLoS 1 ) 


U B 

V b 


W B J 


(3.51) 


(3.52) 


b v 


0(H3(Co) Lg La S T ) 


U 0 " 

V e 

W B _ 


(3,53) 


b V 
6 Z; 


-©(P^N^S 7 ) 


U B 

V 8 

W 0 


(3. 54) 



V 


1 


i-3 

f -i 


a v 
a« 


6(P t L!N T IqS t )- ©(P T N T L t LaSf) 


■1 


V B 

_Wb_ 


(3.55) 


~ ^-©(P’R^Oleiy^R^ + Aoi^s 1 ) 


Ub 

V b 

Lwd 


(3.56) 


«-$(P t N t L,L 3S t ) 


U b 

V 8 

W b 


dv 

a© 


« 0(P t N t La u S T ) 


Ub 

V e 

W B 


(3.57) 


(3.58) 


b v 


= -Q(P t N t LaS T ) . 


U a 

V B 

W b 


(3. 59) 


7* 

If the alternate expression for pa is used (equation (3. 26)). then the largest 


! 

'n 

3 
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tenn is the $ dependent term, so that 
v =a -i (R!f-e 0 )T T La) 


U B 

V b 

Wb 


(3.60) 


In this case, the parttals of v are 


3v 


* -<RiHo)T' La) 


U a 

V B 

W b 


(3.61) 


a v 
Sfo 


:» 4»(L; Ri (-€ 0 )T T Lg) 


Ub 

Vb 

WbJ 


(3.62) 


«» i(RjHc) Ls T t Ls) 


Ub 

Vb! 

Wb] 


(3. 63) 


a v 
ae 


» -4(Rieeo)%H)LiR|(9RaH)l«) 


U B 

V b 

Wb| 


(3.61) 


—• 4(RiHo)T t LoI«) 


U D 

V B 

W b 


(3.65) 


i 
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As in the case of the partials of fD], the "variable" parameters in the 
above expressions for the partials of r v] should be replaced by another set 
of "non-varying" parameters. Then, equation (3.31) can be evaluated for 
the partials of the time delay with respect to the parameters (denoted hy 
vector x). 

The time delay equation, if the small p^, — dependent correction term 
is neglected, and the error in time synchronisation is added, is: 

T ftj *“ (Do • So) + ^tc (3. 67) 

c 


where 6 tc is the clock synchronization error. The Do vector can be 
expressed as 



where 


d is the length of the D 0 vector 

or L , 6 t are the angular coordinates of the vector D 0 
in the inertial coordinate system. 







The changing direction Do in space due to the rotation of the earth changes 
the angular coordinate or L . The vector s 0 is given by equation (3. 14). 

Thus equation (3. 67) gives 

T = — Tsin sin 6o + cos 6 t cos 6 0 cos(a L -ao)] + 5t e . (3.69) 

If the transformation parameters involving precession, nutation, 

polar motion and earth’s rotation rate are assumed known, then the number of 

unknowns are three for the intersite vector, two for the source vector in 

space of each observed radio source and one for the error in initial clock 

synchronization. Hence the total number of parameters for n observed 

sources for a baseline is 2n + 4. However, by assuming the knowledge of 

the rotation axis direction and rotation rate, the origin of longitude is still 

not defined. By fixing one more parameter for the definition of longitude 

origin, the number of independent unknowns reduces to 2n+ 3. Any 

number of observations made to one source from a baseline can only resolve 

— * —¥ 

three quantities which are the components of Dfc>* so along the earth's ads of 
rotation plus 6t c , and the coefficients of the sinusoidally changing components 
of D 0 * So along two orthogonal directions in the earth's equatorial plane. In 
order to determine D, at least three observations must be made of at least 
n sources where n is the smallest integer that satisfies the equation 

3n 2 2n + 3 (3. 70) 

which gives n to be 3. 

The absolute accuracy of the velocity of light in vacuum is about 1 ppm. 
which is much less accurate than the measurement of time delay. However, 
whatever value of c is used in the computations merely provides a scale for 
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the whole system of measurement. 

3.32 Moon-Moon VLBI Equations . 

The equations in this section are derived for the case in which there 

are two or more radio antennas located on the moon. These stations, if 

located at sufficiently large distances, provide a network of selenodetic 

control, whose positions can be determined by the VLBI method. The 

observing procedure would remain basically the same as the earth-located 

VLBI stations although slight modifications may be necessary for practical 

purposes. It is assumed that the measured quantities would be time delays 

as in the case of earth-earth VLBI. 

Let VLBI observations be made at two stations A, B on the lunar 

surface, whose position vector from the lunar center in the 3elenodetic 

- ♦ — * 

coordinate system (see Section 2.32) are r A , r s respectively. 

From Source 

✓ 

KL 


Figure 3.7 Moon-Moon VLBI Observations. 

The components of the two vectors, in the xyz selenodetic system are 
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(3.71a) 


X* 


r A cos b A cos t K 

y» 

= 

r» cos b A sin t. A 

z» 


r A sin b A 


’xa~ 


r a cos bg cos -t B 

Yb 

= 

r e cos ba sin -t- B 

. z ®_ 


r 9 sin ba 


(3. 71b) 


If the vector from A to B is r then 

r = r A - r 8 (3. 72a) 

and in matrix notation 


r » - 


'xa 


X A “ 

r r 

= 

Yb 

- 

y« 

- r *_ 


_Zs_ 


- z »_ 


(3. 72b) 


The vector a is a unit vector from the lunar center to the radio source. 
Because of the small diameter of the moon compared to the distances of 
the radio sources, the vector s also represents the unit vector from the 
radio source in the direction of any point on the lunar surface. 

As was pointed out in Section 3.31, the VLBI time delay (r) can be 
calculated approximately by 


r 


1 _ 

c 


(r 0 • so) 


(3. 73) 


”* -4 — > -* 

where r 0 and So are the vectors r, s whose components are expressed in 
the same coordinate system. The most convenient coordinate system for 
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space orientation on the moon is the ecliptic coordinate system. Conse- 
quently, the r and s vectors can be expressed in an inertial coordinate 
system represented by the 1950.0 mean ecliptic system. 

If the 1950. 0 geocentric mean equatorial coordinates of the source 
is given as a 0 , 6 0 . the geocentric mean ecliptic coordinates are obtained 
by a positive R t rotation through an angle equal to the mean obliquity of 
the equator of 1950.0 tc 0 ). Thus the mean ecliptic coordinates Xo, So are 
obtained from 


Hi(fo) "cos 5 0 sin ar 0 
-sin 5 0 


The "lunar monthly parallax of the moon" as defined by Kolaczek Tsoi is 

very small for the radio sources and the effect of this parallax on 

translation of the center of the ecliptic coordinate system from the geocenter 

to the selenocenter can be neglected. Therefore, the selenocentric mean 

ecliptic coordinates of the source is also given by equation (3.71). 

The r vector defined in equation (3.72b) has its components r„, r r , r t 

along the selenodetic Cartesian coordinate system. This coordinate 

system is related to the mean ecliptic system of date by the three Eulerian 

— * 

angles 0, 0 and 4», as defined in Section 2. 34. The r 0 vector, with com- 
ponents along the 1950. 0 mean ecliptic Cartesian coordinate axes is given 
by 



(3.75) 






where 


P is the precession matrix 


and 

T = H,(*)R t (-d)Ife($) . 


In analogy with equation (3. 1G), the exact equation for the time delay T, 
taking into account the motion of one of the two VLB1 stations between the 
interval r expressed by equation (3. 73) is 

T = j- (Vo • So) [l +■ (3.76) 

where 

ra 0 is the velocity of site B in the same coordinate system in 
-* — * 

which ro and s: are defined. 

The vector r 0o is obtained by differentiating with respect to time, the vector 
* 

reo» which is 

£ 0 = R; (fo) P’lt i(-f)T T r a (3.77a) 

or 



"r Bl " 


* xT 

r 8 0 = 

r * 3 




_ re -\ 


z 8 


(3.77b) 


f 


f 

i. 
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Since Co and r a do not vary with time, 

-► 

* r Bo = Ri(<b)[^4) (1 oP T Ri( - <)T t ) - fi» (Ra(C>)I« R a <-$) 
R 3 (Zi>Ri(-dT T ) + z^p’LaR^-OT') - e(P T L l R i (-e)T T ) 

- T T ) - e(P T R 1 (-f)R 3 (“l/')LiR 1 (0)Ra(-^) 

- *{P T Ri(-f)T I L3)] rS. (3.78) 


Thus if 


r 0 = r x i + r 2 j + r 3 k 


r Bo = r 01 i + r 8a j + r 63 k = v x i + v 2 j + v 3 k 


So ~ Sii + s a j + s 3 k 

where 

**> — > 

i, j, k are unit vector along the X,Y,Z axes defind the 1950.0 
mean ecliptic system, 

then equation foivhe time delay (from equation 3.76) is 


r = “ (r^ + rgS a + r 3 s 3 ) |l + . VlSl ., t .- V ^ S . a * VaSa 


(3.79) 


As in the case of the ear.h-based VLBI, any differences between time delays 
observed and predicted using equation (3.79) can be used in an adjustment 
process to correct the assumed values of the parameters which appear in 
the prediction equation. In order to do this, the partials of r with respect 
to the parameters are needed. 
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To form the partials, equation (3.79) is rewritten in matrix notation 
as 


r = -“[Vs) + 
where the matrices r, s and v are 


r = r 0 = 


■? (rTsV 


- T s)] 


r a 

r 9 _ 


(3.80) 


s = So 


Si 

s 2 

S3 


and 



V 1 


1 

Vi* 

<T 

- . 

V = 

*2 

= 

r&2 


V 3 


*03 


The partials of T with respect to r, s and v are given by equations (3. 28), 

(3.29) and (3. 30) respectively if the matrix r is substituted for the matrix D 

in those equations. Similarly, since the elements of Jie matrices r, s and v 

are functions of parameters (represented by the vector x), the partials of r with 

respect to the parameters xare given by equation (3,31), replacing D with r. 

The partials and will be obtained by differentiating express- 

d x d X ox 

ions for r, s and v. 

Since 
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V -COS 5 0 003 tto 

s = ss = Ri(f 0 ) -cos 6 q sina 0 » 
S 3 — sin Q!o 




-COS COS Qo 
-cos * 0 sin Go 
-sin 


(3.81) 


sin cos 

Ri(e 0 ) sin sino( 0 
-cos 6„ 


(3.82) 


cos 6 0 sin a Q 


= Ri (c 0 ) -cos «> 3 «o 


(3.83) 


and - 77 - , can be obtained by substituting these equations and 

°% ° a 0 dc ° 

equation (3. 29) into equation (3. 31). 

Also, from equation (3.75) 




Hence, taking partials: 


LiRi(C 0 )P’ft 1 (-C)T T r T 


(3.84) 


{ 
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1 


= R^loPThf-cyT r r 


( 3 . 85 ) 


•— = ^v(f 0 )R 3 (So)URa(-«i)R3(2i)Ri(-f)T T r y j 

d°i 


( 3 . 86 ) 


= R 1 (C 0 )P t L 3 Ri(-<)T t r y 


( 3 . 87 ) 


= -Ri(c 0 )P T L l R l (-f)T T r y 

O C 


( 3 . 88 ) 


-B. = -R 1 (gP T R 1 (-()UT T 1 r 



( 3 . 89 ) 


= Ri(f 0 )P T R l (-OR 3 (-^)LiRi(9)R 3 (-^) r, 


( 3 . 90 ) 


-ff = -Ri(e 0 )I*R,<-e>T T I. 3 r y 


( 3 . 91 ) 
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r. 
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: 

f 

«. 

t 

r 

» 

! 


/? 

£ 

& 

t. 

r* ’ 
^ ■- 

&* 

i 


Since 


r. 


r— “ 

*b 


*a 

T r 

= 

y * 

- 

y* 

Jz_ 


_ zs _ 


Z* 


~ = -(R^P^-eKT) (3.92) 

& T h 


= R^fJP^C-OTT . (3.93) 

O To 


As in the case of the earth based VLBI, the coordinates of the two stations 
cannot be solved for simultaneously. However, the differences in coordinates 
can be solved for, or alternatively, one station coordinates can be solved if 
the coordinates of the other station are held fixed. 

The partials of time delay T with respect to the parameters that appear 
in the expression for r* 6o can be obtained in the same way as was in the case 
of the earth-earth VLBI. In equation (3.78), the most significant term is 
the s& -dependent term. By neglecting other terms, the equation for v can be 
approximated by 


v ~ -<$[Ri(c 0 )I*Rx(-f)T»L 3 ] 



Partials of equation (3.94) are: 


-Jr ~ -4»[L 1 Ri(f 0 )P r Ri(-0‘I 7 I«3] 


x 8 

ye 

ze. 


(3.94) 


(3.95) 


IS 
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* ^fR 1 (€o>L3P r Ki{-f)'I f L 3 ] 


(3.96) 


x B 

ye 

Ze 


S v 


^ » 4>IRi(€ 0 )R3(^o) Is Rs (-6i)Ife(Zi)Ri(-e)T T L3] 


x B 

y B 

Zb 




^ «*> -«*Ri(e 0 )P T I a R 1 (-c)T r I 9 I 


x B 

y B 


y B 

_Z B _ 


« *[Ri(€ 0 )P t L 1 Ri(-OT%] 


X B 

ye 

Zb 


-|J ** *[Ri(e 0 )i'Rx(-e)l< J T r L 3 J 


x B 

y B 

Z B 


3 V 


7q « -*[R 1 (f o )P , R 1 (-e)R 3 (-^)LiRa(0)Ra(-«t)I^] 


x B 

y B 

z B 


_^v 


w §[Ri(C 0 ) p 'Ri(-f)' I< L 3 L3] 


X B 

y B 

z B 




(3.97) 


(3.98) 


(3.99) 


(3.100) 


(3.101) 


(3.102) 


(3.103) 
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The partials of v, with the exception of and , are partials with 

3 f(j ^ 

respect to parameters varying from instant to instant. Both in this case 
and in the case of the partials of r, the "variable" parameters can be replaced 
by "non-varying" ones by expressing each parameter as a function of time or 
similar "independent" variable. This procedure has been demonstrated 
in Chapter 2 for the variables involved in this section. 

An alternate set of equations for the moon-moon VLBI can be developed 
if we consider the case in which the orientation of a moon-fixed coordinate 
system with respect to an inertial system (such as the 1950.0 mean ecliptic 
system) is obtained directly by numerical integration of the moon’s equations 
of motion. If 9, 0 , <p represent the orientation (Eulerian) angles 
defining the orientation of the selenodetic coordinate system with respect to 
the 1950.0 mean ecliptic system, then equations (3.75) and (3.77b) for r 0 
and r Bo become respectively: 



where the transformation matrix T 1 is 


(3. 104) 


(3.105) 


T T - *fcH>)R,(9)R3(-i). 

-+ 

Also, r 8o can be obtained by differentiating equation (3. 105) for r^ Thus 
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\ 


"V 


V, " 


~*e 


= 

v a 

= -^(L, t t ) + etRgenLiRi^RG^-^iT 1 ^)) 

ye 

j93_ 


_y 3 _ 


Ze 


(3.106) 


Since the term dependent on i is the dominating term in this equation for 

*B<j» 


~v,~ 

v s 

-(iT T La) 

"V 

Yb 

_v 3 _ 


_Za_ 


(3. 107) 


With these expressions for r 0 and r Bo , the time delay can be computed using 
equation (3.79). The number of parameters involved in computing time 
delay using this method decreases considerably, with only the six initial 
conditions for the numerical integration of 4-, ib, 8, i, li and 6 replacing 
the precessional elements, mean obliquity of date and the orientation 
angles of the moon with respect to the mean ecliptic system of date. 

-* r* 

The partials of r 0 and r 8o are now taken as before: 



(3. 108) 


dr 0 

a e 


RgHjJL.R^Ra*-*) 


(3. 109) 
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1 

J 
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-T t L 3 



J 

» 




-T T La 


Xa 

y B 


Zb 


dv 

bit ) 


4(LqT t La) 


Xfl 

y B 

Zfl 


? 

| 

i 

| 


*e 




Yb 

Ze 



i 

( 

i 

t 

i 

i 

[ 

i 

f 

f 

l 

I 

r 

l 

r 


and 


dv 

b$ 


4>(T T La La) 


'xe 

yB 

Zb 


bv 

br B 


= -4> (T T La) 


■s. 

(3.110) i 


(3.111) 


(3.112) 


(3.113) 


(3. 114) 


(3. 115) 


(3.116) 


(3. 1 17) 
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An approximate time delay equation Tor moon-moon VLBI similar to 
equation (3.69) for the earth-earth VLBI can be derived from equation (3.73) 
as: 

T a ~ [sin ft, sin ftj + cos ft cos 0o cos(Xt - Xq)] + 6t c (3. 1 18) 

where 

ft., X L are the ecliptic coordinates of the baseline vector r 0 
r is the length of the baseline. 

The changing direction of the vector in space due to the rotation of the 
moon, changes the angular coordinate A*., but the rate of change is approxi- 
mately 27. 3 times smaller than the earth-earth VLBI baseline vector. 

If the unknowns to be determined in the moon-moon VLBI observations 
are limited to the radio sources' positions, baseline vectors and time 
synchronization error, then for each baseline observing n sources, the 
total number of parameters are 2n + 3 as in the case of the earth-earth 
VLBI observations. Also, equation (3. 70) is valid for the minimum 
number of sources (three) that should be observed in order to be able to 
solve for the sources' positions and baseline parameters. 

For a good determination of the baseline parameters, the three minimum 
observations to each observed source should be spaced in time as to allow for 
large variation in the direction of the source with respect to the baseline. 

The main factor that accounts for the motion of the source relative to the 
baseline is the moon's rotation, the sidereal period of which is approximately 
27.3 mean solar days. Hence, for a strong solution, observations to sources 
should be made over this period. This requirement is one of the envisaged 
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main difficulties of establishing selenodetic control on the moon with 
moon-moon VLBI observations. 

Another technical difficulty in VLBI observations on the moon is the 
use of accurate time standards and frequency generators, necessary in 
any VLBI observation. While it is required in VLBI measurements that 
the generated local oscillator frequency be highly stable, and time accurate 
to at least one part in 10 u , it is also essential that the VLBI equipment 
for use on the moon be as compact and light as possible, for easy trans- 
portation, Synchronization of the clocks used at moon's VLBI stations 
also presents another technical difficulty. 

One major advantage the moon VLBI stations have over the earth VLBI 
stations is the fact that the moon has no atmosphere. Hence, if VLBI 
observations on the moon become a reality, the ever-present problem of 
refraction, which affects almost all types of observations made on the 
earth, will be absent. 

3.33 Earth-Moon VLBI Equations . 

In tills section, VLBI equations will be derived for the case when the 
two ends of the baseline are located on the earth and moon. As noted earlier 
in this chapter, this form of observation is non-existent at present, but it 
is a future possibility. It is assumed that observations will be made to 
natural radio sources by the two end stations in a manner similar to the 
existing earth-earth VLBI observations. The observed quantities are also 
assumed to be relative time delays. The derivation of the equations in this 
section essentially follows the methods used in the last two sections. 
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Figure 3.8 Earth- Moon VLBI Observations. 


From Figure 3. 8, 

P is the VLBI station on earth with geodetic coordinate «5, X, h 
M is the VLBI station on moon with selenodetic coordinates b, t, r 



\ 





zs^T-a ivr?.-' 



O is the geocenter 
C is the selenocenter 

s is the unit vector in the direction from the source to the origin (geocenter) 
p is the vector from O to the earth VLBI station (P) 

■4 

r is the vec:cr from lunar center (C) to the moon VLBI station (M) 

•4 

D is the vr.rtor from the earth station to the moon station 
M is the vector from geocenter to the moon station 

pmd 

-4 

C is the vector from geocenter to the selenocenter. 

The following vector relationships hold: 

M = C r (3.119) 

-4 -4 -4 

D = C + r - p . (3. 120) 


If the wave front of the observed radio source reaches the station M at time t t 
and reaches tl e station P at ta, r seconds later, then the time delay r is: 


r = — [D • s ] 
c 


1 + 


& • 


( 3 . 121 ) 


This equation tor time delay is valid provided all the vectors are expressed 
in the same coordinate system. The coordinate system chosen is the mean 
ecliptic coordinate system of 1950. 0. 

"4 

The vector D is given by equation (3. 120) as 
D = C + r - p . 

_ -4 «4 *4 

Expressing the vectors C, r and pin the 1950.0 mean ecliptic coordinate 
system we have: 
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(3. 122a) 


Co = Ri(eo)C 


"c," 


XcQ 

^2 

= Ri(<o) 

Ycq 



_Zcq_ 


r 0 = R l (c 0 )P T Ri(-f)Ti r 


r a 

- R^ColP'R^-OTi 

X 

y 

fa 


z 


Po = R v (€o)P T N T S T p 


‘o» 

Pa 

= R,(c 0 )P T N T S T 

"u" 

V 

_Pi_ 


w 


(3. 12 2b) 


(3. 123 a) 


(3. 123b) 


(3.124 a) 


(3. 124b) 


where 

c 0 is mean obliquity of equator at 1950,0, 

Xcq, Y cq , Zcq are the geocentric coordinates of the lunar center in 
the mean equatorial system of 1950.0 (obtainable from 
the JPL's DE-69 development ephemeris), 
x, y, z are the Cartesian coordinates of the moon's VLBI 
station in the selenodetic coordinate system. 
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U, V, W are the coordinates of the earth VLBI station in 
the average terrestrial coordinate system, 

P is the precession matrix, 

N is the nutation matrix , 

S is the matrix necessary to transform coordinates in 

the true celestial system into the average terrestrial 
coordinate system, 

T* is the matrix which transforms coordinates in the mean 
ecliptic system of date into the selenodetic coordinate 
system. 


Th - Raf'i'K) Ri(~6a) RoOPk) • 
In matrix notation, the vector matrix fD] is 


[D] = Da 


( 3 . 125 ) 


The vector s, from the source to the origin (geocenter) whose components 
are expressed in the 1950.0 mean ecliptic system (Sq) can be obtained from 
the cataloged position of the sources. If 5 0 , CKo are the cataloged declination 
and right ascension of a source (given in 1950. 0 mean equatorial system). 


so = Rj(fo) s 


(3. 126a) 
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t 


£ 

7- 

U* 

.s 



"»r 


“-C03 6o COS oto 

9a 

= Ri(fo) 

-cos 6o sin Qo 

_S3_ 


-sin 6o 


1* 

For the M 0 vector, from equation (3. 119) 

-» -* -» 

M 0 = C 0 + r 0 

so that 

■«* ** r* 

Mq = Co + r o • 


(3. 126b) 


(3. 127) 


From equation (3.122a) and (3.122b), 

Co = Ri(fo) C, 


Ci" 



c 3 

= Rl(€o> 

Yc, 

_c 3 _ 




(3. 128 a) 


(3. 128b) 


The numerically integrated lunar ephemeris (such as the JPL's LE-16) 
gives both the position vector C as well as the velocity vector C in the mean 

J-* 

equatorial coordinate system of 1950. 0. Alternatively, C can be obtained 
by differentiating the analytical expressions for the position of the moon as 
given by any analytical lunar ephemeris. The other component of M 0 in 
equation (3. 127), i. e. r 0 ,can be obtained by differentiating either equation 
(3. 123a) or (3. 123b) with respect to time. Thus, similar to equation (3. 78), 



Jc 

'&Ll 
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r 0 = = Ritecjrr'otl^P'RU-e)^ - ^(RatQLs RaHi) 

R3<Zi)Ri(-€)T^ + 2 1 (P t LqR 1 (-€)T«) - C( P T L l R l (-C)Ti) 

- tMP T R,(-C)I* Tj) + 6 m{P t Ri(-€)R3H^L 1 R 1 (0 w )R3(-4h)) 

- fc^R.Hyr;,!*) lr . 


(3. 129) 


In this equation for the r 0 , the dominant term is the term dependent on 4>„. 
Therefore, r Q is approximately given by 


r 0 ~ R^foir-^P'Ril-OT^) 1 r. 
Hence, from equations (3.127), (3.128b) and (3.130), 

X 

Sis I *** Ri(fo) I 


(3.130) 


TMo! = 


l a 

M 3 


*cq 

fcq 

Z 


CQ 


-4> b (p t r!(-€)t;4 > 


. (3. 131) 


The time delay can now be calculated from equation (3.121) as 
T = (D lS , + DaSa + DbSa) [l + ^ + (3.132) 


i. e. , in matrix notation 

T = — [(D T s) + ^ (D T s v T s)] 
c c 

where 


(3. 133) 



Vi 


M, 

fv] = 

v a 

= 



v 3 


i 


115 



It has been mentioned in Chapter 2 that the orientation of both the earth 
and the moon with respect to an inertial system (for example, the 1950.0 
mean ecliptic coordinate system) can be numerically integrated (see also 
Chapter 4 and [78] ). In this case, if 

0c, 0t, 4e represent the orientation (Eulerian) angles of the 

earth-fixed coordinate system (UVW) with respect to 
the 1950.0 mean ecliptic system 

and 

8„, 0 h , 4m are the moon's orientation (Eulerian) angles of the 
moon-fixed (xyz) coordinate system with respect to 
the 1950.0 mean ecliptic system, 

then the vectors p 0 and ro are respectively: 


where 



T c = R3(4 t )Ri(-e E )R3<&> 
Tm = 


(3.134) 


(3.135) 




Consequently, 


= Ri (€o) Y e q + Tj y j - T E r V 


(3.136) 


Also, 


ro= -fWIoTj) + d„(BaHi.)LiRi0 H )R3(-4))-4(THLa)] r. (3.137) 


Since the 4> M term is the dominant term. 


Hence, 


r 0 “^HtTijLij) r . 


M a | m Ri(Co) y co - 4 h( t;l 3 ) y 

_MaJ JZc q _Z 


(3.138) 


(3.139) 


The time delay is still computed from equation (3.132) or equation (3.133). 

The number of parameters that go into the computation of t is greatly reduced 
if the numerical integration method is used to obtain the earth’s and moon's 
orientation angles. 

As noted above differences between predicted and observed values of 
the time delay are used to correct the assumed values of parameters used 
in predicting the time delay. The partials of time delay r are found by taking 
the partials of equation (3.133) as in the other VLBI cases (Sections 3. 31 and 
3.32). Thus, the partials of r with respect to D, s and v are given by equations 
(3. 28), (3. 29) and (3. 30) respectively. Also, equation (3. 31) gives the expres- 
sion for the partials of r with respect to the parameters. The partials of 
[D], [s] and r v] with respect to the parameters will now be formed. 
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From equation (3.129) 


V 

i 


v» 

Da 

= Rl(fo) ! Y c q j + P T 

RiH)Ti : y > - N t S t . 

I i 

V 


— 1 
or 
u 

..Hj 

L Z J 

_ w _ 


This equation for the D vector is identical to the equation for X r vector (for 

laser observations) which is given by equation (2.10) of Chapter 2. Hence 

-» •* 
the partial derivatives of D are equivalent to the partial derivatives of X T 

In Chapter 2 and are given by equations (2.29) through (2.44). The expression 

for the So vector as given by equation (3. 12Gb) is also identical to equation 

(3.74) in Section 3. 32. The partial derivatives of s 0 with respect to Co, 6 0 

and Oo are therefore given by equations (3.81), (3.82) and (3.83) respectively. 

Partials for matrix v is obtained from equation (3. 131) as follows: 


3v 

b € 0 





v 3 



&s 


Xcq| 



Ri (<o) 


(3.140) 


(3. 141) 



-Ri(fo)P r R 1 (-c)TiL 0 ) 



z 


(3.142) 


av 

a Co 





(3.143) 
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dy 

d0I 


= *M[Ri(fo)R3<So) Ls R a (-e 1 )R 3 (z 1 )Ri(-e)TjL 3 ] 


x 

y 

z 


(3.144) 


“ ^-i'sfRt^P'loRU-OT^Ls] 


x 

y 

z 


(3.145) 


-gj = ^[RUCoJP'Li Ri (-OTjLj 


x 

y 

z 


(3.146) 


•fj; - ^IRi(€o)P T Bi(-C)IaTjloJ 


x 

y 

z 


(3.147) 


“fj- = -^B[Ri(fo)P r Ri(-c)R a (-0H)L 1 R l (9 H )R 3 (-*H)L 3 '] 


x 

y 

z 


(3.148) 


= ^«[Rl( f o)P T Rl(“C)T((I/3 1*3 ] 


x 

y 


(3.149) 


dy 

”x" 

y 


_ *«IR:(fo)P T Ri(-c)T5lo9 


(3.150) 


US 





X 


The partials of D and v are slightly different if the alternate expressions 
for them are used (equations (3.136) and (3.139)). Those expressions which 
differ from the ones given earlier are: 


-LsTh y 


(3.151) 


#■ = R 3 (-<|>k)L 1 R 1 (9m)R 3 (-<S>m) y 

O0 M 


(3.152) 


= - TmL 3 y 


(3. 153) 


(3. 154) 


... , ,, . 3D 3D 3D . 3D 

with similar expressions for ^ and 


Also, 


at - 


(3.155) 





(3. 156) 


( 3 . 157 ) 


x 

~ = *„crj i«ifl) y (3.158) 

Z 

V • * 

= -fcrriLo). (3.159) 


t As it was in the case of laser observation equations, most of the 

parameters that appear in the derived partials of time delay in their present 
tV form are "variable” parameters, varying from epoch to epoch. However, 

these "variable" parameters can be replaced by another set of parameters 
which do not vary with time. For example, the integrated Eulerian angles 
(both for the earth and the moon) are functions of the initial conditions and a 
t few physical parameters. These Eulerian angles, which vary with time, can 

V therefore by replaced by the initial conditions and the physical parameters 

gJ 

v" 

r- 
i. • 
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in an adjustment program. This principle has been explained in more detail 
in Chapter 2. 

In the earth-moon VLBI observations involving only one station on the 
earth and another on the moon, if the only unknowns are limited to the inter- 
site vector between the earth station and the moon station and the vector 
to the observed sources, then there are 2n+3 unknowns, where n is the 
number of observed sources. As in the other two VLBI cases, three simul- 
taneous observations each on three radio sources would uniquely solve 
the resulting nine unknowns. 


3.4 Summary . 

In the proceeding sections of this chapter, another new observational 
system - the VLBI - which could be used for selenodetic control determination 
has been treated. The basic operating principles of the VLBI were given, and 
equations relating the geometric time delay to geodetic, selenodetic and 
astronomical parameters were derived. These equations were then differ- 
entiated so as to obtain the observation equations for measured time delays. 
Equations developed in this chapter pertain to earth-earth, moon-moon and 
earth-moon baselines, and are valid only for observations of natural radio 
sources. 

In deriving the observation equations, the observed time delays are 
assumed to have been corrected for atmospheric, instrumental and other 
systematic errors which may affect the measured delays. However, the 
influence of the earth’s atmosphere on VLBI measurements is discussed 
briefly in Appendix B. Also, modified equations for earth-earth VLBI with 
artificial radio sources are derived in Appendix A. 
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4. NUMERICAL INTEGRATION OF THE ORIENTATION OF AN EARTH- 
FIXED COORDINATE SYSTEM WITH RESPECT TO AN INERTIAL SYSTEM 


4. 1 Introduction 

The rotation of the earth on its axis and its orbital revolution around the 
barycenter are of fundamental importance in astronomy. While the circles that 
these motions describe on the celestial sphere serve as fundamental reference 
circles, the apparent motion of celestial bodies produced by the earth's motions 
Is the basis of time reckoning. 

The two motions of the earth are complex and irregular. The rotation of the 
earth on its axis is continually disturbed by the gravitational attractions of the 
sun, moon and the nearer planets. Furthermore, variations in this motion 
occur due to the fact that the axis of rotation does not coincide with a principal 
axis of inertia and because of tidal deformations of the earth. Also, mass 
distribution of the earth is not symmetric, and this distribution changes through 
transfers of mass upon and within the earth from the operation of geophysical 
and atmospheric processes. 

The variation of the direction of the earth's instantaneous axis of rotation 
causes a continuous motion of the celestial poles and the equator among the 
stars, thereby causing variations in the coordinates of celestial objects with 
respect to time. These variations can be divided into two parts, namely the 
secular part (precession), and the periodic part (nutation). Also, there is a 
variation of the position of the instantaneous axis within the earth (polar motion), 
causing a continual, but small variation of latitudes and longitudes 
of points on earth. In addition to the variations in the position of the axis of 
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rotation, there are also variations in the angular rate of rotation, which 
affect the measurement and determination of time. 

In geodetic astronomy, one customarily reduces the cataloged position 
of an observed star to its position at the epoch of observation using the tabu- 
lated values of precessioual elements and nutation in longitude and obliquity 
together with the values of the star’s proper motion. The tabulated values of 
nutation in longitude and obliquity are calculated from, the nutation series 
developed by E.W. Woolard in f 94j. 

The present accuracy to which values of precession and nutation elements 
are known is consistent with the ; evet of accuracy expected from present instru- 
mentation and observations. Furthermore, if the constant of precession is 
in error, the error is absorbed into the stars’ proper motion values due to 
the method by . Inch these proper motions are obtained. Thus, as far as pre- 
sent geodetic astronomy is concerned it would not matter one way or the other, 
whether the value of the precession constant is revised. However, if 
new observational methods of position determination become feasible 
(such as the Very Long Baseline Interferometry applied to geodesy) then the 
need of knowing precession and nutation to better accuracy is of prime impor- 
tance. 

The more precise future observations can be useful in the process of 
finding more accurate values of precession and nutation.' By using a mathema- 
tical model in which precession and nutation quantities are parameters in addi- 
tion to other desired parameters, correction to existing coefficients in the 
expressions for precession and series for nutation can be obtained. The large 
number of the parameters involved (109 for nutation and 9 for precession) 
necessitates a arge number of observations over an extended period of time. 

In order to avoid having to solve for such a large number of parameters, one 
can do either of two things. The first alternative is to cut down on the number 
of parameters by fixing some of the coefficients while solving for the remainder. 


124 



I- 




r.v 

k h 


i. 

r. 

i- 

,y 


V- 

F 

ir 




The second alternative is to devise another way In which the effect of precession 
and nutation could be computed without resorting to series development, and 
in such a way as to drastically reduce the number of the parameters involved. 

The motion of the earth around its center of mass can be represented math- 
ematically by the equations for the motion ofarigidbody. While the earth can- 
not be regarded as perfectly rigid, the effects of departure from rigidity could 
be expected to be Included in the solution for the equations of motion of a rigid 
body obtained by fitting real observations to the mathematical model. 

The three, seconu order differential equations of motion of a rigid body 
have always served as a starting point for finding effects of precession and 
nutation. Understandably in the past, analytical solutions to these equations 
have been derived involving large series expressions. The availability of com- 
puters in the last decade have now made another approach to the solution of 
the differential equations of motion possible. The original equations of motion 
can be integrated numerically without resorting to the various techniques neces- 
sary for finding analytic solutions to the equations of motion. The number of 
parameters involved will reduce to the six initial conditions, which could be 
solved for in a more convenient way than the present coefficients of nutation 
series and precession expressions would allow. By using the numerical inte- 
gration technique, the orientation of any earth-fixed coordinate system (for 
example the "Average" Terrestrial Coordinate System - UVW) with respect 
to an inertial system (for example the Mean Ecliptic Coordinate System of 
1950. 0 - XYZ) could be obtained directly. Thus, the coordinates of an observed 
celestial object desired in the true equatorial coordinate system could be ob- 
tained by first transforming its catalogued position to the UVW coordinate 
system, and using polar motion values and the GAST to transform to the true 
equatorial coordinate system. 

The equations of motion of the earth derived in the following subsections 
is generalized, without simplifications other than the assumption of perfect 
rigidity of the body. 
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4. 2 The Equations of Motion of a Rigid Body. 







t. 






The motion of any rigid body under the influence of external forces can be 
represented as the resultant of a translation, with the velocity of the center of 
mass, and a rotation about an axis through the center of mass. The two com- 
ponents of motion are '.iynamically independent of each other. The center of 

\ 

mass moves as though it were a particle of mass, equal to the total mass of the 
rigid body, and upon which all external forces were directly exerted. The rota- 
tional motion, relative to the moving center of mas3, takes place as if the center 
of mass was at rest. This motion, which will be further investigated here, is 
determined by the moments of the external forces about the center of mass. Parts 
of the equations presented below are derived in more detail in text book? of celes- 
tial and analytical mechanics such as f22T, T85] and T79]. 

Let a rigid body be composed of a system of particles m lf m a m,, . . . , 

whose position vectors from a selected origin of a fixed reference system are i 

■4 ^ ^ -♦ 

respectively r It r 3 , ... r, t .... If external forces F lt F a , ... F 1( ..., and 

4 -♦ — ♦ 

internal forces (due to mutual interactions of the particles), R lt Rg, ... R t , ... 
are the two systems of forces acting upon the system of particles, then the 
force equation for the ith particle is given by 

-* -* d 2 1 , 

F, + R, = (4.1) 

For the entire system, the force equation is given by 

L F ! + E R, = Sm,^ 1 . (4.2) 

The second term on the left hand side of equation (4.2) vanishes since forces due 
to mutual interactions occur in pairs of equal and oppositely directed forces. 

The forces acting on the ith particle is converted to force moments by 
multiplying each term in equation (4. 1) by x, which then gives: 
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•* + ■* H d 3 r, 

r s x F, + r, x R, = m, r, * . 


-» d 3 ?, d ■* dr, 

r “ d? = 7t {T '' di > 


equation (4. 3) gives: 


-* :? -e -* d dr,, 

r, x F, + r, x R, = m, — (r, : l ). 


For the entire system, the force-moment equation becomes 


_ -» _ d -* dr,. 

S r, * F, = Em, - (r, x - ) 


by taking E r, x R, to be zero. 



Figure 4. 1. Rigid Body and a Fixed System. 

Considering Figure 4.1 above, XYZ is a fixed reference system with origin at 
O.and G is the center of mass of the rigid body. If v, is the velocity of the 
particle m, referred to O, then 

-* i* . ;* 

v, = p + r,. 

If the origin O coincides with the center of mass G, then 

v = P = — 1 
Vl dt 


and equation (4. 5) becomes 
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E r, * F, = — [Em, r, * v,J 


L = E r, * F, = force moment 

tt 

H = Em,r, * v, = rotational momentum (moment of momentum). 

n 

For any particular rigid body the force moment L is obtained through the 
force function while H can be expressed in terms of the inertial constants of 
the body and its angular velocity. In the section that follows, the principles 
of rigid body motion will be applied to the case of the earth. 

4.3 The Rotation of the Earth About its Center of Mass. 

4.31 The Dynamical Equations of Motion. 

In Figure 4.2 lettheXYZ axes represent an inertial coordinate system centered 
at O (for example the mean eeliptic system of 1950. 0). Also let the UVW axes 
represent a coordinate system, fixed to the earth's body and also centered at 
O (for example, the average terrestrial coordinate system). The Instantaneous 
axis of rotation OP also passes through the earth's center but does not coincide 
with either the W axis, or the Z axis. 

From equations (4. 6) and (4. 7) 

H = Em, r, * v,. (4.8) 

D 

Replacing the summation sign by an integration sign and writing dm for'the 
element of mass, m, f equation (4. 8) could be written as 

H = f (r, * v,) dm 


where 


r, is the position vector of dm from origin O 
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P — Instantaneous Pole 
u) — Earth's Rotational Velocity 
XYZ — Mean Ecliptic <1950. 0) Coordinate System 
UVW — Average Terrestial Coordinate System 
0 ## — Eulerian Angles 

Figure 4.2 The Fundamental Coordinate Systems and the Eulerian Angles. 
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V 


and v t is the velocity of dm referred to the origin. 


Also, v t could be expressed as 


v, = w * r, 


if m represents the angular velocity of the earth. Hence 

H = J [r, * (w * rj)]dm 

N 

or ^ 

H = J [M(r, . r,) - r,(r, • w)3- (4.9) 

M 

If the coordinates of the mass element dm is u,v,w in the UVW coordinate 
system, and l, j, k are unit vectors along the U,V,W axes respectively, then 

?! = ut + vf + wz 

and 

W = (Ait + m»7 + 

Subtituting for r t and <3 in (4. 9) and carrying out the scalar products we obtain: 
H = J [(u*ji + coyj + M*k)(u 3 + V 3 + wVfui + vj + wk)(UMi +vu>v + wui^jldm. 
Expanding the above equation further and simplifying we obtain: 

H = (Am, - Fw, - Ewji + (BWy-BoJK-FWu)f + (Cu>„-Dtq, - Eu}»)^ 


H = Hyi + H v j + H M k 

where A,B,C,D,E,F are the moments of inertial of the earth defined as: 


A = J* (V s w 8 ) dm 

M 

B = J (u a + w 2 ) dm 

a 

C = J (u a + v^ dm 


D = J vwdm 

N 

E = f wudm 
n 

F = j* uvdm. 


Equation (4. 10) can be rewritten, in matrix notation as 

[H] = [M t ] Cw] 


(4.U) 
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where 


-E -D 


(4. 12) 


In order to differentiate H with respect to time, one also has to consider the 
fact that the directions of the unit vectors i, j, k in space do change with respect to time. 


Thus since 


H = Hyi + H,j + H*k, 


dH _ -fdHu di -dH, di -*dH H dk 
dt - l dt + ^dt + j dt + ^dt k dt ““di' 


From the above expression we obtain 


= (1^ -oj m H v + WvHw) r + (H„ + u>*H, - CAjHJj + (H tt + - w,Hu)k (4. 13) 


which can also be written as 


= M|d) + B» H, 


(4.14) 


This is the right hand side of equations (4. 6) and (4. 7). To obtain the left hand 


side, i. e. L, we have 


as earlier defined. 


L = E ( r t * F,) 


The external forces acting on the earth are mainly due to the gravitational 
attractions of the sun and the moon, and in a much lesser way, to the attractions 
of the planets. Woolard in [94] considered the moon and the sun as the only external 
bodies whose gravitational effects cause precession and nutation. Although the effects 
of the other planets are thought to be negligibly small, the equations to be derived 
here will include their effects so as to make the equations general. Besides the 
generality of the equations, however, it io possible that long-term effects of 
the planets on the rotation of the earth around its center of mass are of such 
magnitude that can be detected with future observational systems. 
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If U is force function and V is potential then 

F = dU = -dV 


.. L = dV * r. 

The potential of the earth of mass M at a distant point P is given approximately 
by the MacCullagh's formula [4l] as 

v = ^ + ^(A + B+C-31) (4.15) 

where G is the gravitational constant, A, B, C are the principal moments of 
inertia of the earth with center of mass at O, R is the distance OP and I is the moment 
inertia of the earth about OP. Equation (4. 15) neglects terms depending on 
1/K* and on, since these quantities are negligible in 
astronomical problems because of 
the great distances involved. 

If we now consider, as 
in the case of the earth, a 
system of n bodies of finite 
dimensions, the total potential 
energy of the system is given by 



where 

Ri* • • Ro are distances from earth's mass center to the 

center of mass of external bodies 1, . . . n respectively 
are masses of external bodies 1, ...n respectively 
A„ Bj, C„ I, are the inertia quantities that pertain to the body i, 

i = 1- • • n In the same way as A, B, C,I were defined above, 
M is the total mass of the earth 
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and I 1 


Is the moment of Inertia of the earth about the line join- 
ing the earth's center to the center of body i, i = 1, . . . n. 




Since V depends on the orientation of the earth only through I 1 , only 
the last term of (4. 16) need be considered. Thus differentiating (4. 16) 


dV = 


3G 


w- 


However, 


I 1 = ^ DclM.x,] 


(4. 18) 


where Xj gives the coordinates of the body i from the earth's center in the UVW 
coordinate system. 

M ( has been defined previously. 

Hence 

r* 

r‘ = 


and 


dl‘ = — a (M, X,] 


dV = -[■~^(M i X 1 1+ 5 |r 3 (M I 3g+ •••+ ^'[M.xjj (4. 19a) 


1 

1 


■t. 


or in vector notation 


Since 


then 


dV = [dV. + dV a + ••• + dV B ]. 


L = dV * r. 


L = dV! * Xi + dV a * X 3 + • • • + dV n « X B 


(4. 19b) 


L = mTXjxXj + ^m7x 2 *X 3 + ... + ^»mTx d * Xb J. (4.20) 


From equations (4. 7) and (4. 14), 

dH 


L = — = M| ta> + w * H . 
dt 
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Therefore, since H = M t u) 



§ 

. x 

V* 



1 -. 

i;. 


i *-• 






-» r3GMs..-*.. 3GM„ 

M, go + aJ * M, (0 = -rr M,X i»X, +••• + - rjr 
- Ri Ro 


M.X »X, 


J' 


(4.21) 


The above equation (4. 21) is the dynamical equation of motion for the earth, 
assuming rigidity, if it i 3 further assumed that the products of inertia (D, E,F) are 
zero, the equation reduces to the well-known Euler's equations of motion [79]. 

4. 32 Euler's Geometric Equations. 

In the last sub-section, the motion of the earth about its center of mass 
has been dynamically represented as a rotation about an axis which constantly 
passes through the center of mass, but whose position within the earth and direc- 
tion in space varies from instant to instant. The orientation of the earth-fixed 
UVW coordinate system with respect to the fixed (inertial) coordinate system 
(XYZ) also varies from instant to instant. This orientation at any instant can 
be defined by a set of three angles - known as the Euler angles - 8, ib, $. 

These angles are presented in figure 4. 2 and are defined as follows: 

8 — the inclination of the average terrestrial equator (moving equator of 
figure) to the plane of the fixed exliptic (of 1950. 0) reckoned positive 
from the plane of the fixed exliptic to the plane of themoving equator, 
such that Q< 90° 

4> — the longitude of descending node of the equator of figure reckoned 
from the fixed equinox 

^ — the angle in the plane of the equator of figure, between the moving 
nodal line and the U-axis, reckoned positive eastwards from the 
descending node to the positive U-axis. 

Thus the relationship between the inertial XYZ coordinate system and the mov- 
ing UVW axes fixed in the earth's body is given by: 


134 






V = Y 


where Rst^), Ri(-6) and Ra(0) are 3*3 orthogonal transformation matrices. 

The rate of change of the three Eulerian angles -4* 0,8 can also be used 

to represent the motion of the earth about its center of mass. In the figure 4. 2, 
• , 

$ is a rotation about the W-axis, 0 a rotation about the Z-axis and 6 a rotation 
about the line of nodes. The relationship between the components of the angular 
veloctiy w along the U,V,W axes — u\j, co y , u) w — and$, 0, § can be obtained, 
by geometrical considerations, from figure 4.2 to be 


<Aj = -0cos$- 0 sin 0 sln<g> 

(t)y = 0 sin $ - 0 sin 0 cos # 

U>„ = 0 cos $ + & 


(4. 23a) 


which could also be written as: 


W, = ^ Hi 


(4. 23b) 


W = SxE 


where 


-cos<P -8ln0 sin© 0 


Sx = 1 sln$ -sinOcos© 0 , E = 0 


Differentiating equation (4. 23c) with respect to time, 


(4.23c) 


• *♦• *)„ 

W = SiE + SxE = Sifi + SaEs 


(4.24) 
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where 


-cos0sin4> sin $ -sinQ cos <t 


cos0cos$ cos sin0sin4> 


-sin© 


E 2 = 84 


From the equations (4, 23) and (4. 24), the relationship between the co and ci> 

r* “> 

vector and the E and E vector is established. These equations, usually referred 
to as Euler's geometric (or kinematic) equations, could now be used in conjunction 
with the dynamical equation (4.21) to obtain three second order differential equations 
of motion, whose solution would give the orientation of the earth-f ixed coordinate 
systenvwith respect to the coordinate system fixed in space,at any particular instant. 

4.33 Differential Equations of Motion of the Earth . 

From the earth's dynamical equations of motion given by (4. 21), expressions 
^ ■> 

for cp and <p given by (4.23) and (4. 24) could be substituted. This gives: 

-* _ -» . r3GM, -> -* 3GM -* -» 

M.SiE + M l S a E s +S v E .MiSiE = - M,X x * X x + • • • + M.X. x X n . 

K n J 


(4.25) 


Retaining only E on the left hand side, 
E = [M.S,] 1 '-3G^ l B mTx x * Xi +• • • +■ 

v LKi 


^IM.X^X^-M.SjEa+M.StE.SiEj . 

(4.26) 


Equation (4.26) is a system of three second-order differential equations of mo- 
tion, and the solution for E gives the motion of the axis of figure - UVW, or 
equivalently, the orientation of the earth in space. The position of the Earth 
relative to the instantaneous axis of rotation at any instant can also be obtained 


1,36 





through the angles that this axis makes with the three body-fixed axes - UVW, 
since the direction cosines of these angles are given by: 




4.4 Numerical Integration of the Earth's Orientation Angles 


The simplified form of equation (4.26) has always been the starting 
equation for the development of the existing solutions of the earth's 
precession and nutation. The differential equations have been derived with 
only one assumption — that the earth is a rigid body. In previous solutions, 
further simplifications are made by assuming that the M (moment of inertia) 
matrix is a diagonal matrix (putting D = E = F = 0). It is also usually 
assumed that the earth possesses rotational symmetry about the W axis 
(i. e., A - B). Furthermore, the external bodies whose attractions influence 
the earth's motion are limited to the moon and the sun, hence the term 
"luni-solar" precession and nutation. 

The simplified differential equations of motion me. prior to this time, 
analytically solved by the well-known method of variation of parameters, 
involving successive approximations. First, the motion that would occur, 
if all external forces were to vanish is solved for, and this solution is there- 
after modified to obtain a solution under actual conditions. The motion of 
the axis of rotation is generally calculated in two parts. The secular part 
of the motion is referred to as "luni-solar precession" and the periodic 
part is called "nutation". A very good example of the classical method of 
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of Woolard [94], which at present serves as the basis for the calculation of 
nutation tabulated in the various astronomical almanacs. 

In this section, an alternate method of using the differential equations 
of motion (equation (4, 26)) to obtain the orientation of the earth with respect 
to an inertial coordinate system will be giver.. The method involves the 
use of digital computers to numerically integrate the earth's equations of 

motion. 

The three second order differential equations (4.26) are transformed 
into a form suitable for numerical integration by setting 


such that 


or 



(4.28) 


where f(E, E) is represented by the righihand side of equation (4, 26). This 
is a system of six first order equations, of which the first three are linear 
and the last three are non-linear. 

These six first order equations can be numerically integrated by any 
standard method. Two methods were used in these studies. The first 
method is the modified Hamming fourth-order predictor-corrector method. 
This is an integration procedure that is considered stable and which uses a 
special Rungc-Kutta procedure to obtain starting values. Provisions for 
controlling local truncation error by changing the step size at each step 
are included. The detailed algorithm is described by Ralston in [80]. 

The second method used is the variable order Adams predictor-corrector 
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method, developed and used by JPL in their Orbital Determination Programs. 
In addition to being a variable order method, provision is made for making 
either absolute or relative error tests for the control of local truncation 
errors for each of the integrated quantities individually. Stepsizes are easily 
halved or doubled at each step. A description of the integration algorithm 
is given by Krogh r 55 1. 

The variable order Adams method was preferred to the Hamming's 
method due to its special features. Although the investigations presented 
in this chapter were started with the Hamming's integration subroutine, 
it was later replaced by the variable order subroutines. Consequently, 
results presented later in this chapter were obtained by using the variable 
order integration subroutines. 

4.41 Calculation of Initial Values . 

Any numerical integration procedure is formulated for initial-value 
problems. It is therefore necessary’ to obtain the starting values: 


Y v, = 


E 1 - 
E lu 


L £ jto 


0o 

Co 

4>o 

00 

Co 

io 


(4.29) 


The starting values were computed in the following way. In Figure 4. 3, 

XYZ — is the mean ecliptic coordinate system of 1950. 0 
UVW — is the "average" terrestrial coordinate system 
0, C> 4- — are the Eulerian angles (orientation angles). 


Then let 




4 

VA 

M 


139 






i, j, k represent unit vectors along XYZ axes respectively 

•4 4 — > 

i , j , k represent unit vectors along the axes UVW respectively. 


and 

n represents the unit vector along the line of nodes, positive 
toward the descending node of the equator. 



Figure 4. 3 Calculation of Eulerian Angles and Their Rates 
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The transformation from the XYZ coordinate system to the UVW 
coordinate system is given by 


~u 


x" 

V 

= SNPRi {-«o) 

Y 

_w_ 

a poch 

Z 


(4.30) 


or 


where 


u 


x” 

V 

= T 

y 

w. 

3 pech 

z 


19S0.0 


T = SNPRi (-Cc) 

and S, N, P and e 0 have been defined previously. Since the 3 * 3 matrix T 
transforms a set of coordinates expressed in the 1950. 0 mean ecliptic system 
to the average terrestial system, it follows that 


T u 

Tia 

Tl3 

Tai 

T aa 

Taa 

T 31 

T 3 a 

T 3 3 


j 

W 


(4. 31) 


From Figure 4 . 3 and equation 4. 31, 0 is given by 
cos 9 = k ' • k = T 33 


Hence 


Similarly, if 


0 = cos l (T 33 ) 


(4. 32) 


= 360° - t}> 
' = ? • n 
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cos 4> J 






k' x k = n sin 8 


n “ '7ITZ C^aai - T ati) 


cos*' = i- " = 


ib' = cos 1 '—£a£ N 
* \sin 8/ * 


(4. 33) 


For the third angle <&, 


cos $ = i ' • n 


cos $ - gin g (An A 3 3 - AiaAox) 


= cos 


(4.34) 


For the initial epoch, the T matrix is calculated. Then the inittal values 
of the three Eulerian angles 8, t and 4> are obtained using equations (4. 32), 
(4.33) and (4.34). 

The time derivatives of the angles are obtained by numerical differenti- 
ation. The angles are calculated for an epoch to + ut and the time rates 
are obtained from 


(4. 35) 


+ M ~ 
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At is a small increment, whose choice is rather arbitrary. However, the 

• > * 

increment must be small enough to make the calculated 9, i j) and § close to 
their true values. From trial runs, it could be concluded that a value of 
1/16384 of a day (approximately 5 seconds) is sufficiently small for this 
purpose. 


4.42 Parameters of Integration . 

In numerically integrating the system of equation (4.28), the function 

for calculating £ has to be evaluated. This function, (equation (4. 26)) is 
-a 

given in terms of E, E and some other parameters. These other parameters 
include the moment of inertia matrix of the earth M )t the masses of the 
moon, sun and the nearer planets as well as their geocentric positions. 

The earth's moment of inertia matrix is given by 


M, = -F 


-E -D 


The elements of the M, matrix are related to the second degree spherical 
harmonic coefficients as follows (see [68]) 

Cao = -KTC - &(A + B)] 

Cai = K • E 


Caa = K 


§ai = K • D 


^2 = K ‘ n 


(4. 35) 


where 


a,, M are the earth's semi-major axis and mass respectively. 




v •* ‘ ■ ' *vr. * v .•* 


• ' # 



From equation (4. 36), there are five quantities which explicitly give the 
six elements of the M, matrix. Consequently, an additional equation must 

be defined in order to find the moments of inertia from the spherical 
harmonic coefficients. From equation (4. 36), 


A - C + ^ (Cso “ 2Cas) 


B = C + — (C a0 + 2C sa ) 


K 821 


(4. 37) 


E " K Cal 


F = ^ Saa 


By putting 


K' = 77ta = C • K 


Ma, 


C = ^ 
^ K 


and equations (4. 37) become 


A = ~ (K' + C so - 2C aa ) 

B = ~ (K' + C 30 + 2C 3a ) 

and expressions for D, E, F remain the same. If by some other means the 

value of K' can be found, then the "scaled" moment of inertia matrix can be 

calculated without using the value of K. It turns out that the scaling of the 

M| matrix by does not affect the validity of equation (4. 26). 

K 
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In the integration program, the value used for K'is taken from 
Kaula T44] 


K' = - ~-5 = 0.33076. 
Maf 


Also, the coefficients of spherical harmonics used were taken from 
Rapp f 81]. Thus, 


M, - - 


0. 32967 -0. 1692 * 10" 

-0. 1692 * Iff 5 0. 32968 

L 0 0 


0 

0 . 

0. 33076 J 


In addition to the sun and the moon, four other planets were used in 
equation (4.26). These are Mercury, Venus. Mars and Jupiter. The 
position of these planets as well as the sun and the moon are obtained in 
the integration program by reading these values from the JPL Development 
Ephemerfs 69 (DE-69) tape. The DE-69 is the latest of JPL's ephemerls 
tapes, and contains the position of all the planets, the moon as well as 
nutation in obliquity and longitude. The positions in the DE-69 were obtained 
by numerical integration process, and the ephemerls is said to be gravita- 
tionally consistent. The planetary masses used are from the JPL system 
of planetary masses F62]. Table 4. 1 gives the gravitational constant for 
each astronomical body. 
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Body 

Gravitational Constants 
(Km7day“ 3 ) 

Earth 

0. 2975542 * 10 16 

Moon 

0.3659906 * 10 14 

Sun 

0.9906936 x 10 sl 

Mercury 

0.1655848 * 10 1B 

Venus 

0.2425068 * 10* 

Mars 

0.3197127 * 10 le 

Jupiter 

0.9458682 x 10 19 


Table 4. 1 

Gravitational Constants 


Adjustment of Initial Conditions 


In order to obtain the earth's orientation angles at epoch t, the six first- 
order differential equations of motton (equation (4. 28)) are numerically 
integrated from a starting epoch Ig. The values obtaim 1 .'or the angles and 
their rates at the desired epoch t depends to a large extent on the starting 
values of these quantities at the initial epoch, and on the parameters that 
are treated in Section 4, 42. The method through which the starting values (at to) 
can be calculated has been given in Section 4. 41 and it necessarily 
involves the current expressions for precession, nutation and Greenwich 
Apparent Sidereal Time, as well as a form of numerical differentiation. 

Small errors in the calculation of the initial values at the starting epoch are 
propagated through the integration from time to to t, and become exaggerated 
over a long period of time. Consequently, it is important to provide a means 






I * • 


by which initial values and parameters can be corrected in an adjustment 
process, involving some form of observations. 

Suppose Z is a set of data which can be computed as a function of the 
earth's orientation angles together with some other parameters, i. e. 


Z = Z(Y, w) 


<4. 37) 


where 


■[51 


and x is a set of parameters also needed in (4. 37). Then the variation in 
Z is obtained through the linearization of equation (4. 37) by a Taylor series 
expansion truncated at the first degree: 


5Z = A T 6 Y + A x 6x 


(4. 38) 


where 


However, the variation 6Y in the Eulerian angles and their rates is further 
related to the variations in the initial conditions Y 0 and other parameters a 
through the solution of the differential equations 


Y - Y (Y 0 ,a) 


(4.39) 


where 


The variation in Y obtained in a manner similar to equation (4. 38),is given 


g|^ 



6Y = U6Y 0 + Vda. 
The matrix of partial derivatives 

~ BE 

U x U 2 " BE 0 

U = 

U 3 U< BE 


BE 

BE 

BE 0 

"iio 

BE 

Bfe 

BE 0 

BE 0 


(4.41) 


is known as the state transition matrix, which describes the transition of a 
differential variation of the initial epoch conditions from time t© to time t. 
Also, the matrix of partials: 


(4.42) 


is called the parameter sensitivity matrix. It describes the effect of a 
differential variation in the parameters a on the integrated quantities E, E. 

When equation (4.40) is substituted in equation (4. 38), the variation in Z 
is obtained as: 


6Z - A v • U • 6 + A, • V • 6a + Ar fix . 


(4.43) 


The matrices of partials A r and A K can be evaluated from the formulas 
obtained by forming the partials of the various parameters which are functions 
of the observables. Such expressions for the laser distances were given in 
Section 2. 4 and for VLBI, inSectlons 3.31, ,3.32 and 3.33. However, the solution of 
the differential equations (equation (4.39)) cannot be written in a closed form. 
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Consequently, approximate expressions for the U and V matrices have 
to be found which can be evaluated by numerical differentiation or inte- 
gration. 

4.51 Evaluation of the State Transition and the Parameter Sensitivity Matrices. 

One way by which the U and V matrices can be evaluated is by numerical 
differentiation (also called numerical variant) method. In this method, the 
desired partial derivatives are obtained by incrementing each of the initial 
conditions, and integrating the differential equations of motion to the 
required time t. Then the nominal values of the angles (obtained using 
the original initial conditions) are subtracted from the values obtained 
when each initial condition is incremented and the result is divided by the 
increment. 

Thus a set of forty-two differential equations of motion can be integrated 
together, six for the nominal initial conditions (Y^) and six for each of 
the six initial values varied. If, for example the first initial condition Y© 
is varied from Y^ to Y^, then 


U=,i 


Y <"fc,> " Y (*o n ) 
" Y ° 


(4.42) 


and similarly for other values varied. 

In these studies, the numerical differentiation method was used to 
obtain the U matrix. The choice of increment is rather arbitrary, although 
it must be such as to give as close approximations to the partials as is 
possible. A number of increments is tried, until by trial and error, the 
U matrix can be considered stable. The stability of the U matrix is ensured 
when a small change in the increment doe3 not affect the U matrix appreciably. 
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The V matrix can be obtained In the same way as described for the 
U matrix. The parameters of interest are varied one by one and the 
equations of motion integrated to the desired time t. Then the nominal 
values of the angles are subtracted from the values obtained by varying 
any of the parameters and the result divided by the increment. 

An alternative way of evaluating U and V is by numerical integration fa]. 
The equations of motion (4.28) can be written in the form 

Y = Y (Y, Of) 

and the variation in Y is given by 

6Y = B5Y + H6a. (4.43) 

Substitute equation (4. 40) into equation (4. 43) to obtain the equation: 

6Y = B* U* 5Y 0 + ('}• V + H) 6a. (4. 44a) 

However, equation (4. 44a) can also be obtained by differentiating equation 
(4. 40) with respect to time: 

6Y = U6Y 0 + ^6a. (4.44b) 

It follows from equations (4. 44a) and (4. 44b) that 

U = B-U (4.45) 

and 

V = B'V + H. (4.46) 

Equations (4.45) and (4.46) are differential equations whose integration 
will result in the state transition matrix U and the parameter sensitivity 
matrix V. The solution to equation (4.45) depends on the solution to the 
equations of motion because the matrix B is a function of the integrated 
values of Y. Once expressions for B and H have been derived, the differ- 
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entlal equations (4.45) and (4.46) can be integrated along with the six | 

equations of motion (4. 28). The matrices B and H are given by j 




Expressions for B al , Ba a , and H 3 can be obtained by differentiating the 
second order differential equations represented by equation (4. 26). This 
differentiation was not carried out in this work since another method was 
chosen for the evaluation of the state transition matrix. 

4.52 Numerical Fit of the Numerically Integrated Angles to the Calculated Angles. 

In order to investigate the feasibility of numerically integrating the 
earth's orientation angles, the differential equations of motion (4.26) 
were integrated for a limited period of time. The results are presented 
in Section 4.6 of this work. As would be expected however, the numerically 
integrated angles differed from the calculated angles. Also, the differences 
grew larger as the epoch of integration moved away from the initial epoch. 

The differences in values of the integrated angles and the calculated angles 
are due, at least in part, to the Inaccurate values of the initial conditions. 
Better initial conditions can be obtained through adjustment processes involv- 
ing real observations together with the use of the state transition matrix and 
the parameter sensitivity matrix already treated in Sections 4.5 and 4.51. 
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Without available real observations, it is still possible to obtain 
corrections to the initial values. In this case, the calculated angles could 
serve as pseudo-observations, while the integrated angles are "approxi- 
mate" values. This idea originated from the work of JPL in the develop- 
ment of their numerical ephemeris. The lunar ephemeris obtained by 
numerical integration was fitted to the Brown- Eckert ephemeris in order 
to adjust the initial conditions for their integration [29]. It is important 
to note that such a numerical fit of integrated quantities to their correspond- 
ing values obtained through other means does not distort the numerical 
integration process. It merely seeks to relate the integrated values to the 
values obtained through other proven means. This is clearly permissible 
in this case because classical methods of calculating the earth's orientation 
are known to yield results very close to the actual orientation of the earth. 

Let 



be the earth’s Eulerian angles obtained through precession, nutation and 
GAST. Also, let 



be the corresponding values from numerical integration with initial condi- 
tions: 
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In a perfect situation, it is expected that 




Since E c is a function of Y<>, it follows that 
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Using the principles of least squares, the weighted sum of squares of AE 
is minimized. Thus, using the notations of Uotila [93!, the corrections, X, to 
the approximate values of the initial conditions Yq is given by 


= ~(A T PAr(A T PL) 


(4.51) 


where 


3(AE) 

»<Yp» 


(4.52) 




L = Mle - (Pb 


<&c - 


(4. 53) 


and P is a weight matrix. In obtaining the design matrix A, the state 
transition matrix is also needed. 


d(AE) _ o(AE) qE 
A " S(Yc) bE ' d(Y 0 ) * 


(4.54) 


However, 


= I (identity matrix) 


■ u 2 ]. 

6(Yo) |_dEo dE u J 1 2J 


Hence the A matrix can be evaluated along with the integration of Y, if the 
Ui, U a submatrices are obtained through one of the methods in Section 4.51. 
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The I matrix is easily obtained by calculating the Eulerian angles at each 
epoch of comparison (using precession, etc,) and subtractingthese values from the 
integrated angles. In general, the size of the A matrix is 3n x 6 and that 
of L is 3n x 1 where n is the number of discrete epochs over which the 
adjustment is to be made. 

The choice of the function to be minimized in a case like this must be 
made carefully. If the quantities are large, it may be necessary to 
find other functions to minimize. The su of squares of the £E quantities 
have been minimized in these studies because they are quite small. Numeri- 
cal results from the adjustment and the improvement it made to the inte- 
gration made with the unadjusted initial conditions are presented in 
Section 4. 6. 


4.6 Numerical Experiments and Results . 

In the previous sections of this chapter, a method by which the orientation 
of the average terrestrial coordinate system with respect to a fixed celestial 
coordinate system (such as the mean ecliptic system of 1950.0) can be obtained 
by direct numerical integration of the earth's equations of motion has been 
outlined. The earth's dynamical (rotational) equations of motion were derived 
in a more complete form, without assumptions as to the mass distribution 
and dynamic shape of the earth. Since the integrated angles and their time 
rates at any epoch depend on the initial conditions and physical parameters 
used in integrating the equations of motion, an outline of a method through 
which corrections to the initial conditions and the physical parameters can be 
obtained was also given. Such a new theoretical approach to an old problem 
as this needs some confirmation a3 to the validity of the method numerically. 

It is therefore the purpose of this section to provide such needed numerical 
support for the theoretical parts of the chapter. 
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An extensive numerical experiment on this subject is well 
beyond the scope of the present study. Consequently, the experiments 
reported in this section are limited to those necessary to verify the correct- 
ness of the equations derived and the computer programs developed for the 
numerical integration of the Eulerian angles and their time rates. Also, 
the ability to adjust the initial conditions with the adjustment computer 
programs based on the method presented in Section 4.5 was verified. 

An important aid in verifying the computer programs developed in con- 
nection with this study was the Simulated Earth-Moon Environment data 
(hereafter referred to in this section as the simulated data) which was 
developed by Papo T 781. The generated simulated data are based on a 
moderately complex model of the earth-moon dyn-mic system consisting 
of a rotationally symetric rigid earth and a perfectly rigid moon whose 
dynamical shape is that of a triaxial ellipsoid. Details of the model and 
the mathematical formulation of the equations of motion of this simplified 
earth-moon system arc contained in T78]. 

A simulated ephemeris of the earth and the moon was created for a 
period of one year beginning at 2440222. 5 JD (1969. 0). The ephemeris 
consists of numerically integrated geocentric position and velocity of the 
moon, the Eulerian angles and their time rates for the earth and the moon 
recorded at half daily intervals. In addition, a fifth-order modified Everett 
interpolation formula f 77 i was available for use in interpolating the 18 
quantities at epochs which fall between tabulated values. 

The following two numerical experiments are reported in this section. 

(1) Fitting the numerically integrated earth's Eulerian angles to 
those obtained from the simulated data. 

(2) Comparing the numerically integrated Eulerian angles to their 
counterparts obtained through classical method (using 




i &. .i—s, w, 
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precession, nutation and polar motion a3 outlined in 
Section 4.41). 
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4.61 Fitting the Numerically Integrated Eulerian Angles to the Simulated 
A ngles. 

As previously mentioned, the availability of the simulated data enables 
one to independently check the equations of motion of the earth as developed 
in this chapter and the computer programs written on the basis of the 
equations. Therefore, the first task that was performed was to compare the 
numerically integrated angles to the simulated angles. In order to make the 
two sets of angles compatible, the numerical integration program 
written 'or the real case had to be slightly modified to accomodate the 
following restrictions imposed by the simulated environment model: 

(1) For the elements of the moment of inertial matrix, A = B and 
D = E = F = 0. 

(2) Only the moon is the external celestial body whose potential 
affects the rotation of the earth. 

(3) The geocentric position of the moon at any epoch is that given 
by the simulated data rather than the one obtainable from a real 
lunar ephemeris. 

It turned out that these modifications to the real integration program were 
easily achieveable through the alteration of a few statements in the program 
and data cards. 

The integrating subroutine used to integrate the equations of motion is 
the BVBQ subroutine, a variable step, variable order Adams integrator f55], 

A comparison of angles integrated and their counterparts from the simulated 
data shows residuals which were less than O.'OOOl for 8 and 0?001 for $ 
and § over the one year interval. These results Indicate perfect agreement 
between the two numerical Integration programs which are results of equations 
derived independently, using different methods and expressed in different forms. 
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The next step was to test the adjustment program and the ability of the 
adjustment method to recover the initial conditions. For this purpose, the 
theoretical initial conditions (values of the integrated quantities at the initial 
epoch of 2440222. 5 JD), read from the simulated data were varied. Three 
test runs were made in which 

(1) only the initial angles were varied 

(2) both the initial angles and their time rates were varied 

(3) only the time rates of the angles were varied. 

The adjustment was performed over an interval of forty days beginning at 
2440222. 5 JD and ending at 2440262. 5 JD. Values of the angles integrated 
with wrong initial conditions were compared with the "true" values at half 
daily Intervals. Table 4. 2 shows the "correct" initial conditions, and the 
initial conditions used in each of the three cases mentioned above. 


Symbol 

Correct 
Initial Values 

Initial Values 
for Case 1 

Initial Values 
for Case 2 

Initial Values 
for Case 3 

0 

0. 4091596226 

0. 4091693189 

0. 4091644707 

0. 4091596226 

<!> 

_-o. mastiff 8 

-0.110152* Iff 4 

-0. 616711 *10 8 

-0. 131898* 15 s 

* 

4. 8950936587 

4.8951033498 

4. 8950985068 

4. 8950936587 

e 

0. 659363 kIO 7 

0. 659363 *1(P 

0. 725299 *10 7 

0. 725299 *1<P 

* 

-0. 300932 *10® 

-0.30093 ‘10 s 

-0. 303942*15® 

-0.303942* Iff 1 

4> 

6. 3003883741 

6.3003883741 

6. 3003943826 

6. 3003943826 


Table 4.2 Starting Initial Conditions for Three Test Cases . 
(9, $ in radians, 0, &, 4> in rad/day) 
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Symbol 

Correct 
Initial Values 

Recovered 
Initial Values 
for Case 1 

Recovered 
Initial Values 
for Case 2 

Recovered 
Initial Values 
for Case 3 

8 

0. 4091596226 

0.4091596227 

0.4091596227 

0. 4091596226 

0 

-0. 131898 xlO 6 

-0. 131869 xl8 B 

-0. 131880 xl0 B 

-0. 131860 xlO 5 


4. 8950936587 

4. 8950936584 

4. 8950936585 

4. 8950936584 

0 

0. 659363 *10 7 

0. 652045 xlO 7 

0. 654780 xlO 7 

0.649609 *10 7 

0 

-0. 300932 xlif 

-0.298465 x10 s 

-0. 299432 x15 s 

-0.300307 xiO 3 

4> 

6. 3003883741 

6. 3003883718 

6. 3003883727 

6. 3003883735 


Table 4. 3 Recovered Initial Conditions for Three Test Cases after Two Iterations . 

(8, 0, <3? in radians, 0, 0, 4> in rad/day) 

Table 4. 3 presents the adjusted values of the initial conditions for the 
three cases after two iterations. Comparison of the recovered parameters with 
their correct values shows agreement up to the ninth decimal figure (equivalent to 
O."OO02 for the angles, and 0."0002/day for their time rates). The residuals 
obtained from the comparison of the integrated angles and those from the simu- 
lated data were, for all three angles and in all three cases, less than O.'OOOl 
throughout the adjustment interval. Since Case 2 is the most general case, the 
residuals before and after adjustments are presented in Figure 4. 4 for 0 and 0 
and in Figure 4. 5 for $. 

From the adjustment results, it appears as if some of the parameters in 
the solution have poor separation from other parameters. This can be seen In 
Table 4.4 which presents the correlation matrix obtained from the adjustment. 

The correlation matrix for all three test cases is identical up to the first two 
significant digits. Also, the inverted normal matrix (weight coefficient matrix) 
for all three cases is very similar, and Table 4. 5 shows the wetght coefficient 
matrix for Case 2. In spite of the high correlation between certain parameters 
(especially between 0 and 4) the results of these experiments show that the 
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solutions for all the parameters do actually converge to their true values. 
The adjustment performed over a 40-day interval with two Iterations takes 
approximately five minutes to run on an IBM 360/75 computer. 



in the Simulated Case. 
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4. 62 Comparison of Numerically Integrated Eulerian Angles to those 
Obtained through Classical Method . 

The experiments performed in Section 4. 61 have demonstrated that the 
method and equations developed in the previous sections of this chapter can 
be used to accurately integrate earth's Eulerian angles and, if necessary, to 
adjust initial conditions in a simulated model. Nevertheless, the ultimate 
test of this method as to its applicability to the real world rests on how 
such numerically integrated angles compare with their counterparts using 
existing methods. The purpose of this section is to evaluate, through 
numerical experiments, whether or not the equations and methods presented 
earlier can be used in a real situation. In this case, there is no known 
"true" solution, and only good estimates of the solution can be obtained, 
usually through an adjustment process. Consequently, the experiments 
reported in this section show: 

(1) How close the values of the angles integrated are, compared to 
corresponding angles obtained through classical methods 

{?.) Results obtained, if the initial conditions are adjusted. 

In Section 4. 41, a method has been outlined through which the three 
Eulerian angles 0, tb, $ can be obtained at any epoch (using values of pro- 
cessional and nutational elements, and the GAST). The initial values of the 
angles for the numerical integration as well as the angles compared with the 
integrated angles at selected epochs were computed using this method. The 
time rates of the initial angles were obtained by dividing the difference between 
computed angles at the initial epoch and at 1/65536 day later by the time 
interval. 

The integration was performed over an interval of approximately one 
year between 2437610.5 JD and 2437970. 5 JD. The earth's equations of 
motion were integrated by the Dv DQ subroutine. The geocentric position 
of the sun, the moon and the planets (see Section 4. 42) at the epoch of 
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Integration were read off the JPL's latest numerical development 
ephemeris — DE69, The constants of integration were those given in 
Section 4. 42. The integration over the one-year period took approximately 
five minutes on the IBM 360/75 computer. 

Instead of the Cowell type of integration performed in the simulated 
case, the Encke type of integration was used in this case. This can be 
done by modifying the original equations of motion from Cowell's to Encke's 
type. i. e. defining a reference case of motion for which analytical expres- 
sions can be given, and obtaining the difference between Cowell's equations 
of motion and the equations of motion of the reference case. Thus, the 
equations integrated are "perturbations" of the reference case. In addition 
to this modification, it was decided to integrate perturbations of the quantities 
9, 0, 0 + 4> and their time rates instead of 9, 0 . 4> and their time rates. The 
main reason for doing this was to find out whether or not the strong correla- 
tion between 0 and 4? can be removed. As it turns out, this affected, very 
little, the high correlation between the two parameters. 

The quantities integrated (in the Encke mode) were: 


where 


6 ; 


"9- 

6 a 


0 - & 

63 

= 

<S>+ 0 - K 0 - K t t 



6 

6 5 


0 

6 s 


4>+ 0 - K, 


9- = 0. 409170 rad. 

= 6 . 280350 rad. 

K 0 = <S»c + 0t = 7. 082407 

K v = 6.30038810 rad/day (sidereal rate) 
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t = time in JD from the initial epoch. 

Also, the quantities compared at each epoch with the computed values were 
0, $ and <t> + tl>. The initial values of the integrated quantities were 

fi! = 7. 00084121 Xl0" 7 rad. 6 4 = -1. 85309594 xl0 _7 rad/day 
6 a = 4. 42734527 xl0" 7 rad. 6 S = -3. 32660736 xl0 _7 rad/day 
6 3 = 1. 30133957 X10 _<i rad. 6 0 = -6. 40369792 x 10 _7 rad/day. 

All quantities integrated refer to the mean ecliptic system of 1950. 0. 

The residuals in 0, ip and $ + $ are presented in Figures 4. 6 and 4. 7. 

It can be seen from these figures that residuals in 0 and $ exhibit a semi- 
monthly period while residuals in <J> +i /> have a dominant seculartrend. However, 
the residuals were very small, not exceeding 0."02 in 8, 0?05 in ij> and 0.'3 In 
0 +■ il> over the period of one year. 

An attempt was also made at adjusting the initial conditions by fitting 
the numerically integrated angles over the computed ones. For this purpose, 
a relatively short interval of forty days beginning at 2437610. 5 JD was 
chosen. The initial values of the integrated quantities (5i to 6 b) were those 
given above. After a single iteration, the corrections to the starting values 
of the initial conditions were less than OfOl and the adjusted values of these 
quantities were: 

= 7. 23518822X 10 -7 6* = -1. 96953572 Xl0" 7 

6, = 4. 99809019 X 10” 7 6s = -3. 21460406 X10" 7 

63 = 1. 30475169* 10" d 6e = -6. 36575208 X 10' 7 

This single solution took approximately three minutes on the IBM 360/75 
computer. 

The correlation matrix as shown by Table 4, 6 shows in general less 
correlation between parameters compared to the simulated case (Table 4. 4). 


165 


















However, the high correlation between 6 b and 6 8 is almost the same as 
that between i and 4 in the simulated case. Table 4. 7 also presents the 
weight coefficient matrix for the adjusted initial conditions. 



Table 4. 6 

Correlation Matrix for Adjusted Initial Conditions 
in the Real Case. 



. Table 4. 7 

Weight Coefficient Matrix for Adjusted Initial Conditions 
in the Real Case. 
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The residuals in 0, 0 and 9+0 before and after the adjustment are shown 
In Figures 4. 8 and 4. 9. It can be seen from these two figures that the 
basic form of the residuals remained after adjustment, but the residuals 
were more or less centered around the zero line. The semi-monthly period 
present before the adjustment remained after the adjustment. 

One basic conclusion that could be made on the basis of tests performed 
In this section is that it is possible to obtain, by the numerical Integration 
method, the orientation of the average terrestrial coordinate system with 
respect to a coordinate system fixed in space (such as the mean ecliptic 
system of 1950. 0). It has been demonstrated that over a period of one year, 
the Integrated angles compare well with the angles computed by classical 
means to within of 2. Also, it is possible to adjust the initial conditions 
through the method given In this chapter in order to obtain a better fit of the 
integrated angles with computed ones using the classical method. Neverthe- 
less, it is important to note that the most desireable goal is to adjust the 
initial conditions, using real observations. 
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5 ADJUSTMENT MODEL FOR LASER AND VLB! OBSERVATIONS, 


5.1 Introduction . 

In the previous chapters, we have treated two basic types of new 
observations which can be used to obtain more accurate position determination 
on the moon. In both Chapters 2 and 3, the prediction equations for the observables 
have been developed as functions of certain physical and astronomical parameters. 
Also partial derivatives essential to the formation of observation equationshave 
been derived. In general, the observed values of the observables and their com- 
puted values will be different, the differences resulting from errors in the 
observations as well as errors in the assumed values of parameters used in the 
prediction equations. These assumed values of parameters can then be 
"corrected" using the deviations cf the observed values from their predicted 
counterparts through an adjustment procedure. Much mo re cften than not, 
there are more equations available than unknowns (due to large number of 
observations), and a unique solution is usually obtained through the use of 
the principles of least squares. This technique Involves the minimization of 
a quadratic form, which is the total weighted sum of squares of the residuals. 

In this study, the two observation types we are concerned with are the 
laser distances and the VLB I time delays. Furthermore, VLBI observations 
have been classified under three types namely 

earth- earth VLBI observations 
earth- moon VLBI observations 
moon-moon VLBI observations. 
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Although the earth- earth VLB1 observations cannot be used directly to 
obtain coordinates of points on the moon (except in the case of observations 
made to artificial radio sources located on the lunar surface), they can be 
used to determine more accurately the orientation of the earth-fixed axes 
in space. In turn, the knowledge of the earth’s orientation in space is 
essential to the determination of selenodetic control through earth-based 
observations. Therefore, the adjustment equations developed in this chapter 
will be applicable to the four types of observation as follows: 

(i) lasers 
(it) earth-earth VLBI 
(111) earth-moon VLBI 
(lv) moon-moon VLBI. 

In developing the adjustment model, consideration is given to the fact 
that values for most of the parameters to be solved for are available through 
some sort of previous determinations. Therefore, since these ’’approximate" 
values of parameters have some variances associated with them, they could 
also be considered as "observations" with corresponding residuals and weights. 
It is also recognized that laser observations are currently under way, while 
the other observation types are only a future possibility. Therefore, a 
sequential solution of the normal equations is desirable, whereby the laser 
observations could be adjusted first, and the result later combined directly 
with other observation types whenever they are available. Lastly, if is 
recognized that there are many parameters involved in the adjustment of this 
kind, some of whlob are either undesirable (nuisance parameters) or known 
more accurately through other means than can be obtained using the observation 
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types considered in this study. For example, parameters related to 
refraction models and clock synchronization errors are nuisance parameters 
while geocentric Cartesian coordinates of earth stations may be obtained more 
accurately through other means, such as artificial satellite data analysis. 

In summary, the adjustment model to be developed will be a generalized 
model in which parameters will be regarded as "observations" associated 
with their corresponding weights. The solution of the normal equations will 
be of a sequential nature, and parameters whose solutions are not required 
will be eliminated from the normal equations while their contribution to the 
solution of the other desired parameters will be accounted for. The derivation 
of the normal equations and their solution follow, in most part, the method 
used by Uotila in [93]. 

5.2 Classification of Parameters in the Solution . 

It has been noted in Section 5. 1 that all parameters in the general ad- 
justment model will be regarded as observations with corresponding weights. 
In addition to these quasi-observations, we also have observations made with 
either laser or VLBI instrumentation. These observations are divided into 
four groups as follows: 

(a) laser — (F - group) 

<b) earth-earth VLBI — (G - group) 

(c) earth-moon VLBI ~ (H - group) 

(d) moon-moon VLBI — (J - group). 

The parameters (or quasi-observatlonc) are also in two categories. 

The first class of parameters are of direct interest, and their solution is 
the main goal of these studies. Rsraraeters in the second class can be 
described as transient. They are needed for the mathematical modeling of 
the observation groups, but they are either of no direct interest (nuisance 
parameters) or they are not as sensitive to the observation types being con- 
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j sidered here as other known types of observations. 

I 

! In this work, parameters are classified into four groups - X 1( 

| Xa, X* - through the consideration of the above defined classes of desired 

[ and non-desired parameters as well as other factors. 
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Parameters Xi 

This group includes the Cartesian (x, y, z) coordinates of the lunar 
stations in the selenodetic coordinate system. These are the main set of 
parameters in which we are interested. Transformation into polar coordinates 
can be done easily (since the lunar figure approximates a sphere). 

Also in the X| parameter group are the lunar orientation parameters. 

In this chapter, the orientation of the moon with respect to an inertial system 
will be assumed to be defined by the three Eulerian angles obtained through the 
numerical integration of the moon's equations of motion. Consequently, the 
orientation parameters consist of the initial values of the Eulerian angles and 
their time rates at a standard epoch. This set of parameters which is non- 
variant, is related to the moan's orientation angle at each epoch through 
a 6 * 6 partial derivatives matrix known as the state transition matrix. The 
pnrtlals of the mathematical functions with respect to this set of parameters 
are obtained by using the state transition matrix together with the partials of 
the functions with respect to the orientation angles at each epoch of 
observation. Thus, the moon's orientation parameters to be solved remain 
six, irrespective of the number of epochs at which observations are made. 

The last set of parameters in this group relate to the dynamic figure 
of the moon. These are the moon's principal moments of inertia, or the 
lower degree and order harmonic coefficients of the moon. These constants 
appear explicitly in the moon's equations of motion which are integrated in 









order to obtain the moon's orientation angles (see [78]). 

All the three sets of parameters in this group are desirable parameters. 
In addition to this, they also appear in the mathematical functions for the 
three observation types F, H and J. 

Parameters Xa 

The parameters in this group are the geocentric Cartesian coordinates 
of the observing stations on the earth and the earth's orientation with respect 
to an inertial system. They feature in the mathematical functions for F, G 
and H-type observations, since these observations require at least one station 
on the earth. 

The orientation parameters are defined in the same way as the moon's 
orientation parameters above. In Chapter 4, the numerical Integration of the 
earth's orientation (Eulerlan) angles and their time rates have been suggested, 
ha this way, the number of parameters reduce to the six initial conditions 
which are related to the Eulerlan angles and their time rates at any observing 
epoch by the 6 * 6 state transition matrix. 

Also included in this group are the parameters of polar motion (coordi- 
nates of the true (instantaneous) pole from the average pole). For easy 
parametisation, it Is reasonable to regard these coordinates as constant over 
short Intervals of time. 

The X a group of parameters may or may not be regarded as desirable 
depending on the purpose of adjustment. In this work we shall treat them as 
desirable parameters. However, i< is clear that the coordinates of stations 
determined from the types of obssn’ations considered in these studies will be 
less accurate than those obtainable with current methods of station position 
determination on the earth. 
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Parameters 

This group of parameters consists of the undesirable (nuisance) param- 
eters which are necessary for modeling observation types F, G and H only. 

The group includes parameters necessary for modeling the systematic errors 
in the existing lunar ephemerls. Also, included in this group of parameters 
are the corrections to the error models used in computing corrections to the 
observed quantities, for systematic observational errors such as atmospheric 
refraction and earth tides (see Appendix B). The parameters will not be 
solved explicitly, but their contributions will be included in the solution for 
the desired parameters. 

Parameters Xa 

This group of parameters is also undesirable, but does appear in the 
mathematical functions for the G, H and J observation types (VLBf). The param 
eters include the radio source positions and clock synchronization errors. Also 
residual instrumental errors (which are systematic in nature) can be included 
in this group. The group of parameters can be treated in the same way as 
the Xg parameters. 

Summary 

In summary, the following tables present the classification of observations 
in the general adjustment model. 


Type (Symbol) 

Description 

F 

Earth to Moon laser observations 

G 

Earth to Earth VLB I observations 

H 

Earth to Moon VLBI observations 

J 

Moon to Moon VLBI observations. 


Types of Observations 




Symbol 


X, 


Parameter Types 

<0 

Cartesian coordinates of 
moon stations 

(If) 

Orientation parameters of 
the moon at Initial epoch 

(Hi) 

Principal moments of 
inertia of the moon 

<t> 

Cartesian coordinates of 
the earth station 

<ii) 

Earth's orientation 
parameters at an initial 
epoch 

(Hi) 

Polar motion parameters 

(i) 

Lunar ephemeris system- 
atic error model 

(H) 

Residual refraction errors 

(Hi) 

Residual earth-tidal errors 

<0 

Radio source positions 

(H) 

Clock offset errors 

(lit) 

Residual instrumental error 


Comments 


Desired parameters 


Desired parameters 


Undesired parameters 


Undesired parameters 


uations. 


The mathematical structure for the adjustment can be written, using 
Uotila's notations (93) as: 


'T »u i 4v-U jUvu. 















F(L?, Lxj» L‘ v L* 3 ) = 0 (5.1) 

G(L£* L 1 ^, L* a t L* 4 ) = 0 (5.2) 

H(L5, L a v L B V L^, L^) « 0 (5.3) 

J(U, L* lt L* 4 ) = 0 (5.4) 

where 


F,G,H, and J are mathematical functions of observations 
ft is a set of laser observations 
Lq is a set of Earth-EarthVLBI observations 
Lh is a set of Earth- Moon VLBI observations. 

Lj is a set of Moon-Moon VLBI observations 

L Xl » L* a , L« 3 , L» 4 are "observation" groups defined in Section 5.2. 

The superscripts denote the following: 

, a - adjusted value 
b = observed value 
o = estimated value. 

Linearization of equations (5.1) to (5.4) through Taylor series yields the 
following corresponding equations: 

ftV, + ft x,V* t + ft * a V* a + ft XgVxa + W, = 0 (5.5) 

+ ft %V, a + ft » 3 V» 3 + ft x 4 V x< + W c = 0 (5.6) 

ft.% + ft H V h + Bh * a V„ a + ft x 3 Vx 3 + Bm ^V* 4 + W« = 0 (5. 7) 

BjVj+BjXjVx, +B J x 4 V» 4 +W j = 0 (5.8) 

where 

Vij, V, a , Vx 3 , V x 4 are corrections to the estimated values of parameters 
V„ V s . V H and V, are observation residuals 
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Br = 4f . Br x, = r|£ W, = F(L,\ L») 


and similar expressions for B^, Bh, B a x,, B„ \V C and \V H . In addition, let P r , 
P s , P H> PjandP Xi represent weight matrices corresponding to L,, 1^., LjUndL^, 
assuming that these groups of observations are not correlated to one another. 

The determination of the corrections (V's) through least squares 
require the minimization of the quadratic form V T PV. Introducing a 
function T and Lagrange multipliers Kr, Kc and K«, it is desired to minimize: 

r = (v/fvv, 4 y/p s v 5 + vJphVh + b (vl p . v. )> 

l~ 1 X 1 1 


-2Kf (BfVr + B,, a V Ia + Bn 3 V, 3 + W ) 

-aK^B, V c + B^V.,* B 5 , 3 V i , 4 ^ 1( V, <+ W t ) 

-2K*’<BhVh 4 BHx t V Xl 4 4 4 b,« 4 V, 4 +W-) 

-2K r J (B J V J 4 B J , l V J , i 4 Bjj 4 Vj 4 ). (5.9) 

Differentiating equation (5.9) with respect to \£ , V c , V„ , Vj and V^, setting 
the result to zero and taking equations (5.5) to (5.8) into account, the system 
of equation (5.10) is obtained in matrix form. This is a system of normal 
equations of the form 

NX + U = 0 

and the solution 

X = -N‘ l U 

usually requires the inversion of the N (normal) matri> 
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5.32 Solution of Normal Equations . 

The solution of equation (5.10) using matrix algebra can be obtained; 
however the dimensions of the normal matrix are large and it is desirable 
to eliminate some of the parameters from the 3ystem. From the first three 
sets of equations, expressions can be obtained for the observation residuals 
Vr, V 6 and Vh as; 

V, = F?Bf Kr 

V 8 = PgBX (5.11) 

Vh = Ph B* Kh 

v, = Pj bJk,. 

By inserting the expressions for V's above into the last three sets of equations 
of equation (5.10) and re-arranging, the system of equations reduces to: 
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Mf 

0 

0 

0 

Brq 


BfXg 

0 


o' 

0 

0 

^*2 

BgXg 

0 

0 

M„ 

0 

Bm Xi 

Bn*g 


0 

0 

0 

Mj 

B *l 

0 

0 

Bf*, 

0 

B l h 

Bjx. 

- p q 

0 

0 



fl5* 

0 

0 

-Pxe 

0 

Bf 


b5% 

0 

0 

0 

_P «3 

0 


Bm H 

BT ^ 

0 

0 

0 


0 ( | V 


Bp Pr’Bf 


B^Bl 

b„p2b; 


= Bj^'b;. 


The equation set (5.12) can be re-written as: 


M B, K 

B\ -P« V, 


where 


0 0 Mj 


Mf o 0 0 Bp Xl Br, 3 Bf^ 0 

0 Mg 0 0 

0 0 Mm 0 ’ * Bhi Bh|j Bmtj Bh^ 


B JXl 0 




Prom the first equation of (5. 13), 

K = -M l (B,V * + W) 
and from the second equation, 

V„ = -(BlM*B* + Pj’ l BlM -l W. (5.14) 

Then, 

L* = Lx + V x . 

Equation (5. 14) above gives the combined solution for the residuals 
of all parameters ("observations'') using all the four types of observations. 
The weight coefficient matrix of L“(adjusted values of parameters ("observa- 
tions")) is given by: 

Qi“x = tP* + BlM'Bxr. (5.15) 

It was already noted above that there are two classes of parameters, 
the first class containing parameters we wish to solve explicitly, while 
the solution to the other class of parameters is not needed. Therefore 
let the V x matrix be divided into two groups namely: 
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V, B - desired parameters 
V, u - undesired parameters 

whore, according to the parameter classification of Section 5.2, 


1 

y 


■p., .i 

L V «.J 

; p « 0 - 

.« ■*%) 

IV 


r»v »i 

1 .< 
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i 

• - 

i 

o 

V 
. ** 

C 


Furthermore, since laser observations will most probably be available 
before the VLBI observations, the solution of V, o using only F-type 
(laser) observations can first be found followed by the sequential solution 
of V, 0 with the addition of the other observation types. 

For the first case (using only F-type observations) the normal equations 
can be written as 



(5.16) 


Eliminating both K, and V„ u from the system, the expression for V» B is 
V* 0 = -[P, B + + Mfi u )^3'« p T , B?« B (Rlf + Mf» u )' l WV (5.17) 

where 


M,« u = 

Also, 

V lg = P, l u ^.K, (5.18) 
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K, = -[Bfr+M^r^v^ + w,]. 


The weight coefficient matrix for the adjusted parameters - L\ 0 Is 
Qi% = p s" P S^^ ,iVV ‘o)'‘ rv -o P - * 

(5. 20) 


where 


M, s = IHP^V 


The last part of the above equation is the contribution to the weight coef- 
ficient matrix, of the ”folded-in” (unsolved) parameter set 

Introducing the G, H and J -types of observations (VLBI), the normal 
equations can be written as: 


M, 0 ^ K, 

0 Me Ke 

^ S -P^ 0 v «* 

0 \ 


Mo 0 0 


Me - 0 Mm 0 ; Be Up - 

0 0 Mj Bj. 


(5.21) 


W. = W M , Kt = | K„ . 
_wj |_K,_ 


The solution of V,, from equation (5.21) is 


\ = ♦M^)*W, ♦(^(S* 


E^)Kkl. 

(5.22) 
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where 







Comparison of equations (5.17) and (5.22) shows that the latter equation 
could be written as 

6V S 


where 


and 


is given by equation (5. 17) 


The term 6 gives the contribution of the G-. II* and J-sets of observations 
if these sets are sequentially used to find after an initial adjustment with 
the F-set of observations. 

An alternative way of carrying out this sequential adjustment is to use 
the system of equations (5. 16) where the F-set matrices are replaced by 
the C-set. However, the "observed" 1^ should In this case be the adjusted 

obtained from the adjustment of the F-set of observations. Also, the 
weight matrix l\ should be changed accordingly to the inverse of the variance- 
covariance matrix resulting from the adjustment of the F-set of observations. 


5. 4 Numerical Exnerlmctas. 

In this study, a few numerical experiments were performed for the 
purpose of investigating the expected accuracies of position determinations, 
using lunar laser distances as observations. The adjustment system is 
based on the laser observation equations derived in Chapter 2 , and the 
adjustment equations presented in the previous sections of this chapter. 
Since no real data were available, laser distances were simulated for 
the purpose of these experiments. No strong conclusions can be drawn 
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from the results of adjustments performed with simulated data concerning 
the results one would ot&ain using real data. Nevertheless, in the absence 
of real observations, simulated data are often used in evaluating a pro- 
posed model, as well as investigating the relative accuracies of the 
various parameters in the model. 

In an earlier report [25], eight lunar control points, arbitrarily, but 
talformly located over the lunar near-side were used to represent positions 
of lunar retro reflectors, to which eight laser stations, located on earth, 
observe lunar ranges. The design of the numerical experiments to be 
reported in this section reflects the existing configuration of lunar ranging 
stations and rctrorcllectors. The three rrtro reflectors deposited on the 
moon by V. S. astronauts serve as the moon stations to distances which 
are measured from laser stations on the earth. A set of numerical 
experiments were performed for the case when the only actively observing 
laser station - the McDonald observatory at Ft. Davis, lexas - ranges 
to the three lunar reflectors. Similar experiments were also performed 
for the case when three laser stations observe laser ranges to the lunar 
reflectors in order to determine any improvement in expected accuracies 
of parameters in the adjustment model as a result of observations from 
multiple laser stations on earth. For this purpose, two laser stations 
were chosen, in addition to the McDonald Observatory, such that the 
stations are located as distant as possible from the McDonald Observatory. 
The two stations arc located at the Crimean Observatory, USSR and the 
Mt. Stromlo Observatory in Canberra, Australia. The establishment of 
permanent laser observing station at these two observatories has also 
been proposed. Tables 5.3 and 5.4 present the coordinates of the three 
laser stations and the three lunar rctrorcflectors. The location of the 
lunar retroOcctors is also shown in Figure 5.1. 

For the limited objectives of the experiments performed, only the 
following parameters were Included In the adjustment model as unknowns. 


i 

3 

$ 
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Point No. 

U(km) 

V(km> 

\V(km) 

Location 

1 

-1330.81462 

-5328.78935 

3235. 69752 

Ft. Davis, Texas 

2 

-4466.54586 

2683.24104 

-3667.44266 

Mt. aromlo, Australia 

3 

3784.28602 

2552.21344 

4440.46175 

Crimean, USSR 


Tabic 5.3 

The Geocentric Cartesian Coordinates of Three Laser Stations 


Point No. 

x(km> 

y(km) 

r.(km) 

Remarks 

1 

1591.42945 

691.97322 

19.26380 

Apollo 11 Mission 

2 

1652. 80823 

-520.44969 

-110.54128 

Apollo 14 Mission 

3 

1554.85951 

99. 26753 

762.35841 

Apollo 15 Mission 


Table 5.4 

Selenodetic Cartesian Coordinates of Three Lunar Retroreflectors 






Figure 5. 1 

location of the Three Lunar Rctroreflectors 

(1> The geocentric Cartesian <U, V, W> coordinates of the 
laser stations on the earth. 

(2> The parameters of the orientation of the U, V, W coor- 
dinate system with respect to the mean ecliptic coordinate 
system of a standard epoch, which was chosen to be 
1969.0 <2440222.5 JO). These parameters are represented 
by the three Euierian angles and their time rates for 
the epoch of 1969.0. 

(3) The sclenodeUc Cartesian <x,y,z) coordinates of the 
lunar retro reflectors, to which laser distances were 
simulated. 

(4) The parameters of the orie&aiion of the x,y,z coordinate 
system with respect to the mean ecliptic system of 1969.0, 
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represented by the physical libration parameters - T, a, o - 

and their time derivatives at the epoch of 1969.0. 

(5) Three physical parameters of the moon, represented by 

the spherical harmonic coefficients C 33 , C a and the 

C-A 

moment of inertia ratio o = — • 

The simulated data were assumed to be completely free from systematic 
errors, and the geocentric coordinates of the moon's center of mass were 
assumed to be known. 

The mathematical model used to simulate the laser distances was that 
presented in Section 2.5, whereby the topocentric coordinates of the lunar 
retroreDector is expressed as (see equation 2.77}: 


(5. 23) 




V 


x* 


V 

1 

j Y„ r 

= 

Yc 

+ P« 

I'm 

- Pt 

V 

|_z„_ 


_Zc 


_Zm_ 


w 


where 

Xc, Y c , Zr are the geocentric coordinates of the center of mass of 
the moon In the 1969. 0 mean ecliptic coordinate system. 
The orientation matrices P t , P* are obtained from integrated Eulerian -ngles 
of the earth and the moon respectively. 

The, two cases of experiments reported in this section are as follows: 

(1) Obtaining the weight coefficient matrix of the adjusted values 
of parameters, using simulated laser distances obtained from 
the Simulated Earth-Moon Environment Data. 

(2) Investigating how well the parameters in the adjustment model 
can be recovered if their "true" values arc altered. This 
investigation was also performed with the simulated environ- 
ment data 
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Simulation of laser distances. 


The flow of computation for the simulation of laser distances are as follows: 

(a) The rotational equations of motion of the earth are numerically 
integrated so as to find the earth's Eulerian angles and their 
time rates at regular half day intervals, over a one-year period. 
All the integrated quantities refer to the initial epoch of 

1969. 0. An additional output of this step is the 6 * 6 state 
transition matrix. The state transition matrix, the Eulerian 
angles and their time rates, and the epoch of computation are 
then stored on a temporary storage device for later use. 

(b) Step (a) is repeated for the moon. 

(c) The U,V,W coordinates of the laser stations, and the x.y.z 
coordinates of the moon points are read as input into a program 
which generates laser distances. 

(d) At each half-day interval, all the "observable" laser distances 
are computed, using information from steps (a), (b), (c) and 

a lunar ephemeris. These distances are then punched out on 
cards for later screening. 

A laser distance is deemed "observable" if the observing 
conditions are such that 

(i) The altitude of the moon point is between 30° and 70° 
(it) The observed moon point is on the front side of the 
moon's disc (with respect to the laser station) and is 
at least 10° off the lunar limb. 

Integration of the Eulerian angles . 

The numerical integration of the earth's Eulerian angles follow the method 
outlined in Chapter 4. The Encke-type integration was used, in which the 
quantities integrated are given as 


3 

ft 

■ l 
s ? 

!i 

n 
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V. 


'b - ik 


* + i/j - K 0 - K,t 


(5.24) 


1>J [_**’+ $ - K, J 

where 0 C , i/>c. K 0 and K ; are constants. The expression for the state transi- 
tion matrix U, which was obtained by numerical differentiation is given as 


(5.25) 


where 


6 = r 6 ( 6s 63 64 6s 63"^' 


6 0 is the value of 6 at the initial epoch of 1969. 0. 

For the moon's Eulerian angles, the equations developed by Papo r 78 1, for 
the integration of the physical libration quantities (r, p, r, ft and o) and the 
state transition matrix (U*), were used. The Eulerian angles are then obtained 
from the physical libration quantities and the expressions for the Cassini laws: 

r*| I"l + v - rf] fr - a] 


(5.26) 









Apart from the physical libration quantities and the state transition matrix, 
the integration program also computes tiie 3x6 parameter sensitivity 
matrix, S, which is given by: 

o = 5 r TU 0 TO p1 T 

QfCaaSCaoV 


A similar expression for the U* matrix is 

u* = a r _T,p p.i.g iZ 

5 r T c CbPc f 0 Ob PoT ' 

•\ 

Weighting. 

The adjustment method used regards all the parameters, in addition to the 
simulated distances, as observations with associated weights. The relative 
weight of each of the parameters was computed as the inverse of its estimated 
variance, thereby choosing the variance of unit weight to be 1. 

The values of the variances of parameters were selected to conform with 
the level of uncertainties in the present knowledge of these quantities. The 
standard deviation of the parameters are as follows: 


U.V.W coordinates, 

m 

= 

25 meters 

x.y.z coordinates. 

m 

= 

1 kilometer 

5j, 5^. 5 .t 

m 

= 

1.0 sec 

6:,. 5 4 , 5 S 

IT* 

= 

0. 5 sec /day 

t, o. a 

m 

= 

20. 0 sec 

f, d. d 

m 

= 

10. 0 sec/day 

Cas 

m 

= 

0. 5 sec 

9 

m 

= 

2 . 0 sec 

Ca-3 

m 

= 

0. 01 sec . 


The standard deviation of each "observed" distance was assumed ' 
to be 15 cm, which is the precision currently expected from lunar 
laser ranges. 


A 
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The purpose of this experiment is to investigate the internal precision 
of the adjustment system. An idea about the internal precision of the 
adjustment system can be obtained through the inspection of the variance- 
covariance matrix, particularly the diagonal elements (variance ), and 
the associated matrix of correlation coefficients. 

The simulated distances used in this experiment were computed using 
the Simulated Environment Data (78}. From laser distances simulated 
according to the procedure already outlined above, a number of them are 
chosen over a certain period of time, for adjustment. There were three 
cases of experimental computer runs: 

(1) In the first case, laser distances were simulated between 

one laser station (at the McDonald Observatory), and the 
three lunar retroreflectors. Different numbers of obser- 
vations over differing periods of time were chosen for 
the formation of the normal matrix. These were as 
follows: 50 observations over a three-month period, 

100 observations over a six-month period, 150 obser- 
vations over a nine-month period, and 200 observations 
over one year. 

(2) In order to have an idea about how the parameter vari- 
ances improve with increase in the period of time over 
which observations are made, the second case consists 
of computer runs where the number of observations was 
held fixed, while the period of observations was varied. 

(3) In the third case, the simulated distances used, com- 
prised of distances between all the three laser stations 
already mentioned above and the three lunar retro- 
reflectors. The total numbers of observations and the 
observation periods chosen were the same as in case (1). 
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The partials of distances with respect to the parameters were evaluated 
using the equations developed in Section 2.5. Since the orientation param- 
eters (of both the earth and the moon) defined in Section 2.5 are slightly 
different for those in this adjustment model, a simple transformation 
had to be made as follows: 


and 


where 


and 


6D 

3D 

3 E C 

o ^ 

d 5 0 

3E C ' 

3 5 

3 5o 

6 D 
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oE m 
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d E* 
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0 ?0 


• - — for the earth's orientation angles 


for the moon's orientation angles 


5 = % 6 a 6 3 5, 5 5 5*V 

E e = r « t tit tit *\ it <n] T 

£ - r T C O T & 6 V 

E — r 6« *i)« 9* * 4 m! 


D is the distance. 

It follows form equations (5.24) and (5.2G) that 
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The expressions for the parlials of distances with respect to the orientation 
parameters can thus be evaluated. 

Table 5.5 gives the diagonal elements of the variance-covariance matrix 
of case (1), where the only laser station observing is located at the McDonald 
Observatory. The general pattern shows, as can be expected, that the param- 
eter variances decrease as more observations are made over longer periods 
of time. It can also be seen from the table that, in addition to a general 
improvement in lunar position determination through the use of laser ranges, 
we can also expect an improvement in the determination of geocentric 
positions of stations on the earth's surface. A general improvement in 
the accuracy of determining the orientation parameters of both the earth 
and the moon, and; the moon's lower order spherical harmonic coefficients 
can also be expected by using lunar laser ranges. The exception to this is 
the variances obtained for the moon's physical libration in the node and its 
time derivative (o and &). Even with 200 observations over a one year 
period, the variance obtained for cy and a are 250'.'7 and 29.5 secVday 2 
respectively, compared to the a priori values of 400" and 100 sec 2 /day 3 . 

A good explanation for this poor determination of cx and a cannot be found, 
especially since the dominant term of libration in node has a period of 
about one month. On the other hand, it should be realized that the node 
itself is defined by the intersection of two planes at an angle of only 1?5. 
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Another trend that could be noticed in Table 5.5 is that, in general, 
the x coordinates of the lunar reflectors are better determined than the 
y and z coordinates. The z-coordinates are the ones with the poorest 
determination. This phenomena can be explained through the fact that the 
ranges are likely to be more sensitive to changes in the X coordinates 
of lunar stations since the X-axis is always oriented approximately toward 
the earth. This same conclusion cannot be made for the geocentric 
coordinates of points on the earth which are likely to depend more on the 
geometry of the simulated distances chosen for adjustment. 

Table 5.6 displays the variances of parameters for case (2), where 
the number of observations are held fixed and the periods of observation 
are varied. All the parameters show sensitivity to the length of the 
period of observations as can be seen from decreases in the variances 
of parameters when the period of observation is increased. The sharpest 
change in the variances occurs when the period of observation is increased 
from three months to six months. A further increase in the period to 
nine months shows further decrease in the variances, but to a lesser 
degree. 

In Table 5.7, the diagonal elements of the variance-covariance matrix 
(variances of parameters) for case (3) are tabulated. In this case, there 
are three observing laser stations on the earth. The characteristics of 
this table are similar to those of Table 5.5 for case (1), where only one 
laser station is ranging to the lunar retroreflectors. From both tables, 
i; can be seen that the accuracy of determining the parameters increase 
with increase in the number of observations and the period of observations. 
The largest degree of improvement is obtained when the number and 
period of observations are increased from 50 and three months to 100 
and six months respectively. Also, the least sensitive parameter in the 
adjustment model remains the moon's physical libration in node and its 
time derivative (o and a). 
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50 Observations 

50 Observations 

50 Observations 

100 Observations 

' \ 

100 Observations 

Parameters 

3 Months 

6 Months 

9 Months 

6 Months 

9 Months 


meter 3 

meter" 

meter 3 

meter 3 

meter 3 

U 

312 

3.03 

1.2 

1.0 

0.43 

V 

19.6 

0.2 

0.1 

0.07 

0.03 

W 

36.0 

30.0 

10.7 

8.9 

2.88 

Xl 

16.3 

5.4 

4.9 

2.9 

2.0 

Vi 

57.1 

16.9 

13.4 

10.1 

7.1 

Zl 

6*5.3 

563.1 

241.7 

220.9 

117.4 

*3 

10.5 

3.7 

3.1 

2.3 

1.6 

y 3 

138.7 

48.0 

35.4 

27.6 

16.8 

z 3 

627.4 

531.2 

211.5 

191.9 

92.2 

*s 


16.3 

14.4 

14.4 

12.3 

ya 

113.2 

45.1 

28.0 

26.7 

19.8 

_ 23 

539.6 

437.2 

168.0 

171.5 

85.5 


second 3 

second 2 

second 3 

second 3 

second 3 


0.6 x 10~ 3 

0.3 x 10" 3 

0.7 x 10"* 

0.1 x 10" 3 

0. 2 x 10"* 

6a 

0. 3 x 10~ 3 

0.5 x 10' 3 

0.7 xlO" 3 

0.2 X 10" 3 

0. 1 x ID" 3 

63 

0.5 

0.4 X10' 2 

0.2 xlO" 3 

0.1 X 10" 3 

0.6 x 10" 3 

r 

1.5 

0.48 

0.4 

0.3 

0.2 

CT 

397.6 

392.4 

372.2 

382.0 

348.5 

P 

1.4 

1.38 

1.2 

1.2 

1.0 

Ca 

0.3 x 10" 3 

0.4 X 10" 3 

0.2 XlO" 3 

0.2 X 10" 3 

0.9 X 10"* 

0 

0.4 x 10" 3 

0.4 X 10"* 

0.3 xlO"* 

0.3 X 10"* 

0. 2 X 10"* 

Ca> 

0. 1 x 10" 3 

0. 1 X 10" 3 

0.1 xlO" 3 

0.1 XlO" 3 

0. 1 x 10" 3 


IESSZ3ISHHM 

in 'Fi i' r i 

mm 

sec 3 /dav 3 

sec" /dav 2 

6* 

0.2 x 10” 3 

0.5 X 10" 3 

0.1 xlO" 3 

0. 2 x 10" 3 

0.2 x 10" 3 

65 

0.6 X 10“ 1 

0.2 x 10' 1 

0.6 x 10” a 

0. 7 X 10" 3 

0. 8 x 10" 3 

«6 

0.4 x 10 -3 

0. 1 x 10" 3 

0.4 x 10"* 

0.5 x 10"* 

0.6 * I®" 6 

T 

0.8 x 10" 3 

0. 1 x 10" 3 

0.5 X 10"* 

0.6 x 10"* 

0.2 x 10"* 

0 

93.2 

94,5 

80.9 

81.3 

68.9 

p 

0.02 

0.015 

0.014 

0.015 

0.013 


Table 5.6 

Variances of Parameters tor Case (2). 
Period of Observation Varied 
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Parameters 


50 Observations 
3 Months 


100 Observations 
6 Months 


Ui 

Vj 

Wi 


meter 3 


169.43 

10.61 

9.69 


moter~ 


1. 17 
0.08 
2.49 


U a 

V 3 

W, 


43.04 

119.04 

8.08 


0.34 
0.81 
2. 14 


u 3 

38. 88 

0. 25 

v 3 

85.47 

0.59 

W 3 

9.83 

2.51 

x, 

10.50 

2.51 

y> 

43.26 

6.92 

Zl 

243.70 


Xa 

6.71 

2.08 

y a 

81.66 

22.11 

Za 

207.30 

96.73 

x. 

17.05 

13.90 

y 3 

67.06 

17.71 

z 3 

197.62 

98.69 


second 3 

second 3 


0. 14 x 10” 3 

0 . 38 x 10~ 4 

6 a 

j 0. 11 x 10" 3 

0.39 xlO” 4 


| 0.25 

0.17 xl0” a 

T 

i 0.94 

0.23 

a 

| 397.58 

378.84 

P 

1 1.44 

1. 14 

Caa 


0.20 xlO" 3 

0 


0.25 XlO -4 

Cgo 

j 0.1 xlO' 3 

0.1 XlO -3 

1 IsecVday 2 

se c 3 /day a 

6* 

■ 0.50 x 1Q~ 3 

0. 18 xlO -3 

65 

j 0.99 x 10" 3 

0.18 xlO* 3 


0.67 x 10"* 

0.12 xlO* 4 

f 

I 0.55 xlO" 3 

0.55 xlO -4 

O 

' 97.37 

| 78.40 

P 

1 0. 16 x 10" 1 

j 0.15 xlO" 1 































The correlation coefficient matrix presetted in Table S.S is that of 
Case (3) with 200 observations over a one year period. The correlation 
matrix obtained for the other cases are similar to the one presented 
here, but, in general, correlation between parameters are higher for 
cases involving less observations performed ever less periods of time. 

The first striking feature of the correlation matrix (Tabic 5.8) is the 
high correlation that exists between the C and V coordinates of the same 
lsser station on the one hand, as well as the high correlation between the 
x and y coordinates of the retroreflectcrs cn the other. Also, the U 
f coordinates are correlated among one another, and the same applies to 

f the U, W, x, y, ?. coordinates. These phenomena may be due to the fact 

f that the orientation of both the earth and the iroon, as well a9 the geo- 

I comer are part of the parameters in the solution. It can also be noticed 

r 

[ in Table 5.8 that there is a high correlation between and 6 a , as well 

| as between and 6 4 . This was also the case in Chapter 4 (see Tables 

l 4.4 and 4.6). For the moon's orientation parameters, high correlation 

* • • 

£ exists between Q and p, and p and O. The coordinates of laser stations 

I 

f were neither correlated with the selenodetlc coordinates of the reflectors, 

| nor with the oriematioa parameters of the moon. 

I Experiment (2). 


t 

h 


r 


This experimem was performed in order to investigate how well the 
parameters in the adjustment model can be recovered, if known shifts are 
introduced into the starting parameters (as approximate values of the 
parameters). The computational methods in this experiment are similar 
to that of Experiment (1), except that arbitrary shifts were introduced 
into the "known" values of parameters. The experiment was performed 
for the case when three stations are observing to the three retroreflectors, 
with only 50 observations over a three month period. 
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Table 5.9 presents the shifts Introduced into the various parameters, 
and the corrections obtained from the solution. The coordinates of the 
laser stations were recovered to within 10 meters in almost all the cases 
while the coordinates of retrorefioctors were recovered to within 70 meters 
for most of them. The earth's Eulerian angles are well recovered 
although their time rates (except t a ) are not. Except for T and p, the 
physical libration parameters are poorly recovered, especially or and b. 
The moon's physical parameters C 22 , 0 and Caj were all well recovered. 

The degree of recovery of parameters for this experiment, especially 
for the coordinates of the laser stations and retroreflectors are some- 
what lower than might be expected. Nevertheless, it is not utterly sur- 
prising because of the high correlation that exists between some of the 
parameters. Moreover, only one adjustment cycle was performed for the 
case presented here, and two or more iterations on this solution might 
have brought the magnitudes of the corrections closer to the respective 
shifts introduced. 
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Shifts Introduced to Parameters, and Corrections from Solution 


206 










G. SUMMARY AND CONCLUSIONS. 


The main purpose of this study has been that of deriving the observation 
equations necessary to utilize the lunar laser ranging and the VLBI measure- 
ments for the establishment of a primary control network on the moon. The 
control network consists of coordinates of moon points in the selenodetic 
Cartesian coordinate system, which is fixed to the lunar body, centered at 
the moon's center of mass and has its axes coincident with the three principal 
axes of the moon. The control points are those moon points to which range 
measurements have been made from terrestrial stations, or those points 
which have been used in conjunction with terrestrial or other lunar stations 
for VLBI observations of natural radio sources. 

The problem of determining coordinates of points on the moon in the above- 
defined selenodetic Cartesian coordinate system using earth-based observa- 
tions is rigorously tied to the knowledge of the following dynamic behavior 
of the earth and the moon: the orbital motion of the moon about the bary- 
center, the rotational motion of the moon on its axis and the motion of the 
earth about its center of mass. In addition, the knowledge of the geocentric 
positions of the terrestrial stations is essential. Therefore, it can be expected 
that our knowledge of the parameters related to the above phenomena of the 
earth-moon dynamic system as well as the geocentric positions of earth 
stations can be improved simultaneously with the determination of coordinates 
of lunar points. For this reason, the derived observation equations for the 
laser ranges and the VLBI time delays were based on a general model, in 
which the following groups of unknown parameters were included: 

1. The selenodetic Cartesian coordinates of lunar points. 



2. The geocentric Cartesian coordinates of earth stations. 

3. The parameters of the orientation of the selenodetic coordinate 
system with respect to a celestial coordinate system fixed in 
space. These may be given by three Eulerian angles, obtainable 
either through the numerical integration of the moon's equations 
of motion, or through the Cassini laws and the physical libration 
angles determined through the classical method. 

4. The parameters of the orientation of the average terrestrial coor- 
dinate system (U,V,W) with respect to a fixed celestial coordinate 
system. Traditionally, these have usually been given in terms 

of the precessional elements, the nutations in obliquity and in longi- 
tude, the Greenwich Apparent Sidereal Time and the motion of the 
true pole with respect to the CIO pole. It has been shown in this 
study that the orientation of an earth-fixed coordinate system with 
respect to a fixed celestial system can be represented by three 
Eulerian angles, which can be obtained through numerical integration 
of three second-order differential equations. 

5. The geocentric coordinates of the center of mass of the moon as 
given by a lunar ephemeris. 

The laser and the VLBI equations were developed in Chapters 2 and 3. In 
Chapter 4, some effort was devoted to the derivation of the earth's equations 
of motion (rotational) and an outline of how these differential equations can be 
numerically integrated was presented. The justification for the extensive 
treatment of this subject in this study lies in the fact that the knowledge of the 
earth's orientation in space is an important factor for accurate determination 
of selenodetic coordinates of points on the moon. Moreover, the number of 
parameters defining the earth's orientation is greatly reduced when the angles 
are obtained through numerical integration process. This eases up the impor- 
tant practical problem of computer storage in handling the solution of a huge 
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system of normal equations, which would have resulted from the inclusion 
of precessional and nutational parameters in the adjustment model. A 
general adjustment model, based on the theory of adjustment computations 
using matrix algebra, was developed for the analysis of laser and VLBI 
observations in Chapter 5. 

A few numerical experiments were performed in this study, basically for 
the purpose of checking the models and equations developed in parts of this 
work. The first group of experiments are limited to the numerical verifi- 
cation of the equations and methods developed for the numerical integration 
of the earth's orientation angles, and the computer programs based on these 
equations. From the results obtained, it was demonstrated that, over a 
period of one year, the numerically integrated angles compare very well to 
the classically computed ones, to within 0"2 in the worst case. It was also 
shown that the solutions for the initial conditions of the integration do con- 
verge inspite of high correlations between some of the parameters. 

Another group of experiments was performed for the purpose of 
investigating the accuracies of determining geodetic and selenodetic 
coordinates as well as the orientation parameters of the earth and the 
moon, when laser ranges are used. For this purpose, the three retro- 
reflectors on the moon deposited by U. S. astronauts served as the 
lunar bench marks and simulated ranges were obtained to their reflec.. ..a 
from three laser stations located on the earth's surface. Results of the 


experiments indicate significant improvement in the accuracy of deter- 
mining station positions (both on the earth and on the moon). Also, the 
results indicate that the orientation of the earth and the moon, represented 
by 6 parameters each, can be more accurately solved for with the use 
of laser ranges. However, it was clear from the results that a good 
number of the parameters in the solution were highly correlated. This 
fact seems to have affected the degree to which parameters could be 
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One logical extension of the study performed here is the numerical verifica- 
tion of the observation equations for the VLBI, and the improvement that might 
be expected in the adjustment system if VLBI observations are combined with 
laser ranges. Another problem area is the investigation of the possibility of 
conbining the types of observations treated in this work with other types such 
as earth-based optical observations, satellite-borne photography, lunar 
laser altimetry and range/range-rate observations to lunar satellites. 

In deriving the observation equations for the laser and the VLBI, as well 
as the equations of motion of the earth, both the earth and the moon were 
assumed to be rigid bodies. The effects of non-rigidity of the earth and 
the moon on these equations still remain to be investigated. 

Finally, the ultimate goal in the application of the findings contained in 
this study to Geodetic Science and other related academic disciplines would 
be reached when the equations presented here are successfully applied to 
actual laser and/or VLBI observations. 
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APPENDIX A 

Equations for Earth-Based VLBI Observations of an Artificial Source 

The use of a radio beacon, deposited on the moon, as an artificial 
source for an independent clock interferometric system was one of the 
recommendations of the Lunar Science and Exploration Conference of 
1967 r75"l. The radio source could be a transponder, whose transmitted 
signals are observed by two or more radio telescopes widely separated on 
the earth. The main advantage of using an artificial radio source is the 
stronger signal reception expected because of the proximity of the artificial 
source. The differences in the arrival times of the radio signal at the two 
interferometer stations (time delay) are obtained through the usual VLBI 
technique of delay mapping (described in Section 3.2). The difference between 
the observed and the predicted time delays can be used to obtain corrections 
to approximate or assumed values of parameters used in the prediction 
equation. 

Since the source is at a finite distance from the earth, the equations 
for the earth-based VLBI observations of natural sources derived in 
Section 3. 31 is no longer valid, and a corrected equation will be derived 
in this appendix. 

In Figure A. 1, 

O is the geocenter 
C is the selenocenter 
A, B are the VLBI stations on earth 
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-» 
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-» 
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-* -» 
R*. Re 
-» -» 
S*. Sfe 
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is the position of the radio source on the moon 
is the vector from geocenter to selenocenter 
is the vector from geocenter to the radio source 
are vectors from geocenter to A and B respectively 
are vectors from A, B to the radio source 
is the vector from A to B . 
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Figure A. 1 Earth-Based VLBI with Artificial Source . 


From the figure, 

MF = MB = }% |. 

Hence, the delay distance (d) is given by AF. If c is the velocity of light, 
the time delay T, is given by 


d 


However, 




d = AG + GM - MB 
or 

d = |D | cos 9 + |Sb { cos y - [Sd | 

and since 

|Sb | sin y = |D | sin 9, 

it follows that 

d = |D | cos 9 - |d| sin9 ( 1 ~^” )« 

Therefore, 

t= ^[cos9- sine( J ^“)] 

represents the corrected equation for the time delay, which is the dif- 
ference in the arrival times (tg, tj) of the radio signal at stations B and 
A respectively. 

If the various vectors in Figure A. 1 are expressed in the same coor- 
dinate system (for example, the mean ecliptic system of a standard epoch), 
then the following vector relationships hold: 

-♦ 4 

D = Rg - R* 

S„ = C + r - R, 

= C + r - Re. 

Also, 

-» -> 

SvJjP 

|s»1|d| 


(A. 2) 


(A. 3) 


(A. 4) 


cos 9 = 


(A. 5) 



** % *-%-r>^.*^ 4 *rTa »v'*-''’?:*rj > »vJr*V-,»-''‘ v1 *'* 


■ V" - 


cosy 


-* -» 

= & - A . 

&l£i 


In matrix form equations (A. 5) and (A. 6) can be written as 


cos 0 = 


(D T DSl^ 


cosy = 


Sg Sa 

(SlS^Sl^ 


^ 4 *♦ 

where D, S Af S& are column matrices of vectors D, S A and respectively. 

The equation for the time delay (equation (A. 4)) can be evaluated for 
any epoch through vectors D, S» and S*. As noted earlier, the components 
of these vectors must be expressed in the same coordinate system, such as 
the mean ecliptic system of any standard epoch. 

The vector D which is the inter-site vector, is a function of the geo- 
centric position vectors of the two observing stations and a transformation 
matrix which varies with time. Both S» and S*, depend on the geocentric 
position of the lunar center of mass, the selenocentrlc position of the radio 
source and transformation matrices that vary with time. 
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APPENDIX B 

Systematic Errors Affecting Earth-Based Lunar Observations 

Most of the systematic errors which have effects upon all observations 
made to points or objects on or near the moon can be classified under the 
following two main groups: 

A. Instrumental errors 

B. Errors due to physical effects. 

In the development of observation equations for the two types of observations 
considered in the main chapters of this work (i.e. the laser and the VLBI). 
it was assumed that these observations have been corrected for all systematic 
errors. This assumption is a reasonable one, in view of the fact that most 
observers go to a great deal of effort to correct their observations for known 
instrumental and other systematic errors. Nevertheless, it often happens 
that inadequacy in the theoretical knowledge of the physical phenomena 
involved and Incomplete identification of all sources of instrumental errors 
may result in Incomplete systematic error modelling and the treatment of 
observations, thereby leaving residual errors In the corrected observations 
presumed to be free from all systematic errors. Consequently, the correct 
procedure in adjusting any set of observations should include provisions for 
the inclusion of mathematical models capable of absorbing all or most of 
the residual systematic errors. While the development of such models for 
laser and VLBI observations is beyond the scope of this present study, 
the purpose of this appendix is to identify and discuss possible sources of 
systematic errors in the types of observations already mentioned above. 
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A. Instrumental Error*. 

Instrumental errors are those errors which are Inherent In the 
measuring systems. Their sources are not easily identifiable in new and 
complex observational systems. Full identification of all the sources of 
errors comes with experience in the usage of such new measuring systems, 
and through the analyses of measured quantities. However, once the sources 
have been identified, instrumental errors can be controlled or corrected very 
accurately. 

Since both the laser and the VLB1 are new observational systems, all their 
sources of errors have not yet been identified. Nevertheless, a few sources 
of errors can be mentioned tor each system. 

For the laser, these sources include errors in the time interval measur- 
r ing circuitry, and in the determination of the physical point in the transmitting/ 

1 receiving telescope, to which the measured distances are referred. Another 

: source of instrumental error in the laser is in the determination of total 

: instrumental delay of the signal, due to the optical path length of the telescope 

| and the delay through the photomultiplier. 

; The main source of instrumental error for the VLBI Is in the possible 

drift of the generated local frequency standard, and the limited phase stability 
of the locally generated signal. This results in a clock offset error which 
j can be easily Incorporated into the mathematical adjustment model as one of 

! the parameters to be solved. Also, inaccurate determination of the "system 

( zero bias" is one of the sources of VLBI instrumental errors. The "system 

i zero bias" is dependent on the determination of the focal point of the telescope, 

• and on the instrumental delays of the observed radio signal. 

Most of the instrumental errors are expected to be eliminated from the 
j observations, leaving only residual errors which, if large enough, can be 

| determined by the incorporation of suitable mathematical models Into the 

s “ 

| 

\ adjustment model. 
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B. Errors Due to Physical Effects. 

Systematic errors In this group are caused by the fact that earth-bound 
observations to moon points are made through the earth's atmosphere, and 
that certain geophysical phenomena (such as earth-tides) do change the 
positions of earth stations relative to the geocenter. Therefore, errors 
in this group can be treated under the following headings: 

(1) Atmospheric refraction errors 

(2) Earth tidal effects 

(3) Continental drift 
1. Atmospheric refraction errors. 

Atmospheric effects on radio or optical signals can be classified under 
the two following headings: 

(a) Tropospheric effects 

(b) Ionospheric effects. 

The lower part of the atmosphere ( the troposphere' contains both dry air and 
water vapor, while the ionosphere (the upper part of the atmosphere) Is made 
up of layers of free, electrically charged particles (ions). The electron 
densities of the ionospheric layers vary’ considerably with solar activity, geo- 
graphy and the time of the day. 

The refractive index in both the troposphere and the ionosphere departs 
from unity, thereby causing both a deviation ot the ray from the straight line 
path, and a change in the phase velocity of the signal from its velocity in 
vacuum. Thus: 



where 

v p is the phase velocity in the atmosphere 
c is the velocity of electromagnetic radiation in vacuum 
and 

n is the refractive index of the medium. 
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The departure of the ray path from a straight line is usually very small for 
normal observing conditions, and it is non-existent in the zenith direction. 
Neglecting this deviation, the difference between the path in free space and 
the actual path is given by: 

5d= I d h — **ndh. 

O v 


or 

6 d = J(l-n) dh (B.2) 

where the integration is to be performed over the entire path. 

The above equation (B. 2) is exact for r.n electromagnetic signal observa- 
tion in the direction of the zenith, and corrections for observations made at 
other elevations are functions of the correction in the zenith direction, pro- 
vided no horizontal gradient is present at the time of observation. 

Tropospheric effects. 

In the troposphere, the refractive index (n) is greater than unity. The 
correct evaluation of equation (B.2) requires the knowledge of the refractive 
indices along the entire ray path at the time of observation. Since this is 
practically impossible, numerous methods have been proposed for the 
prediction of the value of the integral equation. These methods include radio- 
optical dispersion techniques, range measurements at two optical wavelengths, 
and prediction from surface weather data. Hopfield in [35] shows how, with 
the use of surface weather data, the dry air component of the range correction 
for laser and radio observations, in the zenith direction can be predicted « ith 
a root mean square error of few millimeters (out of a total correction of 
about 2. 3 meters). However, the range error prediction for the more vari- 
able but smaller wet part of the troposphere is not as accurate as that for the 
dry part. The residual tropospheric refraction error is larger for the VLBI 
observations, than for laser ranges because of the increased influence of 
water vapor in the radio spectral region. At optical frequencies, the effect of 
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water vapor is so much decreased as to be almost negligible. 

Ionospheric effects . 

Unlike the troposphere, the ionospheric medium is such that its refractive 
index is less than unity. Consequently, the two effects (tropospheric and 
Ionospheric) tend to cancel each other. For a microwave radio path which 
traverses both the troposphere and the ionosphere, the overall effect is a 
negative differential phase path length. 

The calculation of the ionospheric refraction error, which depends on 
the frequency of the wave, the electron density, the earth's magnetic field 
and the collision frequency can be done using ionospheric models and the 
ray -tracing technique f6l]. In general, the ionospheric correction is 
inversely proportional to the square of the signal. From r45l, an approxi- 
mate equation for the average day-time ionospheric correction for the 
observed range is given by: 

|6d| =» cm (B. 3) 

where 

5 d is range error 
and 

f is propagation frequency. 

The observations made at radio frequencies are much more affected by 
the ionosphere than those at optical frequency. Consequently, it is desirable 
that VLBI observations be made at high frequencies, subject to other limiting 
conditions. 

From the above discussion, it is evident that corrections due to atmo- 
spheric effects are usually calculated from theoretical models. Residual 
refraction errors are almost certainly bound to remain in the corrected 
observations., the magnitudes of which depend on the departure of the actual 
atmospheric conditions from the model used in calculating the corrections. 
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VLBI observations involve two stations, and the residual refraction error in 
the corrected time delay depends on the residual refraction errors at both 
stations. The determination of the residual errors, if large enough can be 
done by introducing suitable mathematical models into the adjustment 
model, and regarding these residual errors as unknowns in the adjustment. 

2. Fai th tidal effects . 

The tidal phenomenon, which is a result of the varying attractive forces 
of the moon and the sun, produces periodic deformations of the geoid known 
as earth tides. These deformations cause small variations in the intensity 
of gravity, in the direction of the vertical, and in gcoidal heights. On an 
absolutely rigid crust, the periodic variations in the magnitude and direction 
of the gravity vector can be calculated from the potential of the bodies which 
causes these variations. Theoretically, the magnitudes of the variations 
are Or04 and 0!'2 milligals for the deflection of the vertical and the intensity 
of gravity respectively. However, the globe as a whole cannot be truly regarded 
as rigid, since the earth is somewhat elastic and viscous. Consequently, the 
theoretical deformations of the earth due to luni-solar attractions and their 
phases are altered by the earth's elasticity and viscosity. Furthermore, the 
non-rigidity of the earth gives rise to tensions and cubic expansions which 
would not be present with a rigid globe. 

In lunar laser ranging and earth-based VLBI, the observations are made 
with instruments rigidly attached to the earth's crust. The periodic variations 
of the radial distances of observing stations introduce systematic errors into 
both the laser and the VLBI observations, which should be eiliminated before 
adjusting the observations. 

From the static theory of tides, the theoretical height of the observed 
tide (£) at a point is given by r G3 1 : 

5 = ^ <B.4> 
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\V 2 is the potential (to the second degree) of the force-producing bodies 
(sun and moon) at the point 
g is the value of gravity. 

For the moon, if z is the geocentric altitude, the height of the theoretical 
tide is 

f . = 2G. 7 (cos 2 ► 1/3) cm. (B.5) 

The geoidal height therefore increases by a maximum of 35. Gem and decreases 
by a maximum of 17. 8cm, giving a total range of 53.4 cm. For the sun, 

* = 12. 3 (cos 2 z w * 1/3) cm (B. 6) 

which is equivalent to an oscillation from lG.4cm to -8.2cm or a total range 
of 24. G cm. Thus the combined luni-solar effect would give an oscillation of 
of 78cm for the geoid, if the earth were a perfect liquid. 

From observations of tides, the actual response of the earth to the tidal 
generating forces results in a rise and fall of the elevation of the earth's 
surface on the order of 25 to 30cm in the midlatitudes and of 50cm near the 
equator. The projected precision of both the laser and the VLB1 measurements 
shows that for both instruments, the systematic errors due to earth tides can- 
not be ignored, and observations made bv both instruments must be corrected 
for tidal effects. 

In the classical method of determining tides, earth tides are broadly- 
classified under three types namely: 

(i) long periodic 

(ii) diurnal 

(iii) semi-diurnal. 

Within each type, there are several tidal waves of differing amplitudes and 
frequencies. The instruments used in earth tide observations include hori- 
zontal pendulums, tidal gravimeters and extensometers. From a continuous 
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record of observations made with these instruments, the amplituc j and 
phase of each wave is determined, using methods of harmonic analysis 
which are based on the theory of the conibii.ation of ordinates r G3T The 
laser and the VI.13I can play important roles in the determination of earth- 
tides in future provided measurements made by them are of a continuous 
nature. 

It is relevant to mention here that solid tides on the surface of the moon 
do exist, and that tidal effects on the moon are relatively larger, at least 
theoretically r G3\ So far, very little is known about the moon's internal 
constitution, and no observations have been made on the surface of the moon 
for the purpose of determining lunar tides. Consequently, lunar tides for 
the present, can only be considered theoretically. 

3. Continental drift . 

During the past few years, some geodetic evidence has been accumulated 
which shows that the earth's surface is mobile, l.arge blocks of the earth's 
outer layer appear to be moving with respect to one another at varying rates 
from I cm/year to as high as 15 cni/year r 45l. This "drifting of the conti- 
nents" causes changes in the geocentric positions of earth stations, to which 
laser and VLDI instruments are rigidly fixed. However, the average rate 
of motion is small enough to have no appreciable effects on lunar laser rang- 
ing and VLDI measurements, if the observations cover a relatively short 
span of time. For observations taken over a long period, the effect of 
crustal motion may become appreciable, and the observations should there- 
fore be corrected before they are used for adjustment purposes. 
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